diff --git a/Makefile b/Makefile index 05d90238..831d6302 100644 --- a/Makefile +++ b/Makefile @@ -9,13 +9,13 @@ INCDIRS=-I. -I$(LIBDIR) MODOBJS= psb_prec_type.o psb_prec_mod.o MPFOBJS=psb_dbldaggrmat.o psb_zbldaggrmat.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_dilu_bld.o \ + psb_dmlprc_bld.o psb_dsp_renum.o psb_dbjac_bld.o psb_dilu_bld.o \ psb_dprecbld.o psb_dprecfree.o psb_dprecset.o \ psb_dbaseprc_bld.o psb_ddiagsc_bld.o psb_dgenaggrmap.o \ psb_dprc_aply.o psb_dmlprc_aply.o \ psb_dbaseprc_aply.o psb_dbjac_aply.o\ psb_zasmatbld.o psb_zslu_bld.o psb_zumf_bld.o psb_zilu_fct.o\ - psb_zmlprc_bld.o psb_zsp_renum.o psb_zilu_bld.o \ + psb_zmlprc_bld.o psb_zsp_renum.o psb_zbjac_bld.o psb_zilu_bld.o \ psb_zprecbld.o psb_zprecfree.o psb_zprecset.o \ psb_zbaseprc_bld.o psb_zdiagsc_bld.o psb_zgenaggrmap.o \ psb_zprc_aply.o psb_zmlprc_aply.o \ diff --git a/psb_dasmatbld.f90 b/psb_dasmatbld.f90 index 0e84f36e..6c469abf 100644 --- a/psb_dasmatbld.f90 +++ b/psb_dasmatbld.f90 @@ -54,7 +54,7 @@ Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dasmatbld Implicit None ! .. Array Arguments .. @@ -77,6 +77,7 @@ Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) & tot_recv, ircode, n_row,nhalo, nrow_a,err_act Logical,Parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err + name='psb_dasmatbld' if(psb_get_errstatus().ne.0) return info=0 diff --git a/psb_dbaseprc_aply.f90 b/psb_dbaseprc_aply.f90 index 1cf071f0..d2c06809 100644 --- a/psb_dbaseprc_aply.f90 +++ b/psb_dbaseprc_aply.f90 @@ -41,7 +41,8 @@ subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dbaseprc_aply + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -61,20 +62,6 @@ subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err - interface psb_bjac_aply - subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type), intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbjac_aply - end interface - name='psb_dbaseprc_aply' info = 0 call psb_erractionsave(err_act) diff --git a/psb_dbaseprc_bld.f90 b/psb_dbaseprc_bld.f90 index d4dee7f2..f1d8ad59 100644 --- a/psb_dbaseprc_bld.f90 +++ b/psb_dbaseprc_bld.f90 @@ -37,7 +37,8 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dbaseprc_bld + Implicit None type(psb_dspmat_type), target :: a @@ -46,56 +47,6 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) integer, intent(out) :: info character, intent(in), optional :: upd - interface psb_diagsc_bld - subroutine psb_ddiagsc_bld(a,desc_data,p,upd,info) - use psb_base_mod - 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_ddiagsc_bld - end interface - - interface psb_ilu_bld - subroutine psb_dilu_bld(a,desc_data,p,upd,info) - use psb_base_mod - 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_base_mod - use psb_prec_type - implicit none - - type(psb_dspmat_type), intent(inout) :: 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_base_mod - use psb_prec_type - implicit none - - type(psb_dspmat_type), intent(inout) :: 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,ictxt,& & me,mycol,np,npcol,mglob,lw, mtype, nrg, nzg, err_act @@ -112,7 +63,7 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_baseprc_bld' + name = 'psb_dbaseprc_bld' if (debug) write(0,*) 'Entering baseprc_bld' info = 0 @@ -139,17 +90,10 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) ! ! 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',& & diag_,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) @@ -175,7 +119,7 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) goto 9999 end if - case (bjac_,asm_) + case(bjac_,asm_) call psb_check_def(p%iprcparm(n_ovr_),'overlap',& & 0,is_legal_n_ovr) @@ -188,59 +132,16 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) 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' + if (debug) write(0,*)me, ': Calling PSB_BJAC_BLD' if (debug) call psb_barrier(ictxt) - 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 (debug) call psb_barrier(ictxt) - 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(debug) write(0,*)me,': Done umf_bld ',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 Bjac/ASM??' + call psb_bjac_bld(a,desc_a,p,iupd,info) + if(info /= 0) then info=4010 - ch_err='Inconsistent prec f_none_' - call psb_errpush(info,name,a_err=ch_err) + call psb_errpush(info,name,a_err='psb_bjac_bld') goto 9999 + end if - case default - write(0,*) 'Unknown factor type in baseprc_bld bjac/asm: ',& - &p%iprcparm(f_type_) - info=4010 - ch_err='Unknown f_type_' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end select case default info=4010 ch_err='Unknown p_type_' diff --git a/psb_dbjac_aply.f90 b/psb_dbjac_aply.f90 index 81ab3cc4..96817b03 100644 --- a/psb_dbjac_aply.f90 +++ b/psb_dbjac_aply.f90 @@ -43,7 +43,8 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dbjac_aply + implicit none type(psb_desc_type), intent(in) :: desc_data @@ -63,7 +64,7 @@ subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err - name='psb_bjac_aply' + name='psb_dbjac_aply' info = 0 call psb_erractionsave(err_act) diff --git a/psb_dbjac_bld.f90 b/psb_dbjac_bld.f90 new file mode 100644 index 00000000..b631b90a --- /dev/null +++ b/psb_dbjac_bld.f90 @@ -0,0 +1,368 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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. +!!$ +!!$ +!***************************************************************************** +!* * +!* This is where the action takes place. * +!* ASMATBLD does the setup: building the prec descriptor plus retrieving * +!* matrix rows if needed * +!* * +!* * +!* * +!* * +!* some open code does the renumbering * +!* * +!* * +!* * +!* * +!***************************************************************************** +subroutine psb_dbjac_bld(a,desc_a,p,upd,info) + use psb_base_mod + use psb_prec_mod, mld_protect_name => psb_dbjac_bld + implicit none + ! + ! .. Scalar Arguments .. + integer, intent(out) :: info + ! .. array Arguments .. + type(psb_dspmat_type), intent(in), target :: a + type(psb_dbaseprc_type), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + character, intent(in) :: upd + + ! .. Local Scalars .. + integer :: i, j, jj, k, kk, m + integer :: int_err(5) + character :: trans, unitd + type(psb_dspmat_type) :: blck, atmp + real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8 + logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false. + integer nztota, nztotb, nztmp, nzl, nnr, ir, err_act,& + & n_row, nrow_a,n_col, nhalo, ind, iind, i1,i2,ia + integer :: ictxt,np,me + character(len=20) :: name, ch_err + character(len=5), parameter :: coofmt='COO' + + + if(psb_get_errstatus().ne.0) return + info=0 + name='psb_dbjac_bld' + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt, me, np) + + m = a%m + if (m < 0) then + info = 10 + int_err(1) = 1 + int_err(2) = m + call psb_errpush(info,name,i_err=int_err) + goto 9999 + endif + trans = 'N' + unitd = 'U' + if (p%iprcparm(n_ovr_) < 0) then + info = 11 + int_err(1) = 1 + int_err(2) = p%iprcparm(n_ovr_) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + endif + + call psb_nullify_sp(blck) + call psb_nullify_sp(atmp) + + t1= psb_wtime() + + if(debug) write(0,*)me,': calling psb_asmatbld',p%iprcparm(p_type_),p%iprcparm(n_ovr_) + if (debug) call psb_barrier(ictxt) + call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& + & blck,desc_a,upd,p%desc_data,info,outfmt=coofmt) + + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_asmatbld') + goto 9999 + end if + + t2= psb_wtime() + if (debug) write(0,*)me,': out of psb_asmatbld' + if (debug) call psb_barrier(ictxt) + + + select case(p%iprcparm(iren_)) + + case (1:) + + ! + ! Here we allocate a full copy to hold local A and received BLK + ! Done inside sp_renum. + ! + + call psb_sp_renum(a,desc_a,blck,p,atmp,info) + + if (info/=0) then + call psb_errpush(4010,name,a_err='psb_sp_renum') + goto 9999 + end if + + t3 = psb_wtime() + if (debugprt) then + call psb_barrier(ictxt) + open(40+me) + call psb_csprt(40+me,atmp,head='% Local matrix') + close(40+me) + endif + if (debug) write(0,*) me,' Factoring rows ',& + &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 + + select case(p%iprcparm(f_type_)) + + case(f_ilu_n_,f_ilu_e_) + + call psb_ipcoo2csr(atmp,info) + + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + call psb_ilu_bld(atmp,p%desc_data,p,upd,info) + + if (info/=0) then + call psb_errpush(4010,name,a_err='psb_ilu_bld') + goto 9999 + end if + + + if (debugprt) then + + open(80+me) + + call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') + write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m + do i=1,p%av(l_pr_)%m + write(80+me,*) i,i,p%d(i) + enddo + call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') + + close(80+me) + endif + + + case(f_slu_) + + call psb_ipcoo2csr(atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + call psb_slu_bld(atmp,p%desc_data,p,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='slu_bld') + goto 9999 + end if + + case(f_umf_) + + call psb_ipcoo2csc(atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csc') + goto 9999 + end if + + call psb_umf_bld(atmp,p%desc_data,p,info) + if(debug) write(0,*)me,': Done umf_bld ',info + if (info /= 0) then + call psb_errpush(4010,name,a_err='umf_bld') + goto 9999 + end if + + case(f_none_) + info=4010 + call psb_errpush(info,name,a_err='Inconsistent prec f_none_') + goto 9999 + + case default + info=4010 + call psb_errpush(info,name,a_err='Unknown f_type_') + goto 9999 + end select + + + + call psb_sp_free(atmp,info) + + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + + + case(0) ! No renumbering + + + select case(p%iprcparm(f_type_)) + + case(f_ilu_n_,f_ilu_e_) + + call psb_ipcoo2csr(blck,info,rwshr=.true.) + + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + call psb_ilu_bld(a,desc_a,p,upd,info,blck=blck) + + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_ilu_bld') + goto 9999 + end if + + + if (debugprt) then + + open(80+me) + + call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') + write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m + do i=1,p%av(l_pr_)%m + write(80+me,*) i,i,p%d(i) + enddo + call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') + + close(80+me) + endif + + + case(f_slu_) + + atmp%fida='COO' + call psb_csdp(a,atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_csdp') + goto 9999 + end if + call psb_rwextd(atmp%m+blck%m,atmp,info,blck,rowscale=.true.) + if (info == 0) call psb_ipcoo2csr(atmp,info) + if (info == 0) call psb_slu_bld(atmp,p%desc_data,p,info) + if(info /= 0) then + call psb_errpush(4010,name,a_err='slu_bld') + goto 9999 + end if + + call psb_sp_free(atmp,info) + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + + case(f_umf_) + + atmp%fida='COO' + call psb_csdp(a,atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_csdp') + goto 9999 + end if + call psb_rwextd(atmp%m+blck%m,atmp,info,blck,rowscale=.true.) + if (info == 0) call psb_ipcoo2csc(atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_ipcoo2csc') + goto 9999 + end if + + call psb_umf_bld(atmp,p%desc_data,p,info) + if(debug) write(0,*)me,': Done umf_bld ',info + if (info /= 0) then + call psb_errpush(4010,name,a_err='umf_bld') + goto 9999 + end if + + call psb_sp_free(atmp,info) + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + + case(f_none_) + info=4010 + call psb_errpush(info,name,a_err='Inconsistent prec f_none_') + goto 9999 + + case default + info=4010 + call psb_errpush(info,name,a_err='Unknown f_type_') + goto 9999 + end select + + case default + info=4010 + call psb_errpush(info,name,a_err='Invalid renum_') + goto 9999 + + end select + + t6 = psb_wtime() + + call psb_sp_free(blck,info) + if(info/=0) then + call psb_errpush(4010,name,a_err='psb_sp_free') + goto 9999 + end if + + if (debug) write(0,*) me,'End of ilu_bld' + + call psb_erractionrestore(err_act) + + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + +end subroutine psb_dbjac_bld + + diff --git a/psb_dbldaggrmat.f90 b/psb_dbldaggrmat.f90 index 0801bafd..4dd0d93b 100644 --- a/psb_dbldaggrmat.f90 +++ b/psb_dbldaggrmat.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dbldaggrmat + implicit none type(psb_dspmat_type), intent(in), target :: a @@ -49,6 +50,7 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info) logical, parameter :: aggr_dump=.false. integer ::ictxt,np,me, err_act character(len=20) :: name, ch_err + name='psb_dbldaggrmat' if(psb_get_errstatus().ne.0) return info=0 diff --git a/psb_ddiagsc_bld.f90 b/psb_ddiagsc_bld.f90 index 0e01356d..2b8e0d3c 100644 --- a/psb_ddiagsc_bld.f90 +++ b/psb_ddiagsc_bld.f90 @@ -36,7 +36,7 @@ !!$ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_ddiagsc_bld Implicit None type(psb_dspmat_type), target :: a @@ -61,7 +61,7 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_diagsc_bld' + name = 'psb_ddiagsc_bld' if (debug) write(0,*) 'Entering diagsc_bld' info = 0 diff --git a/psb_dgenaggrmap.f90 b/psb_dgenaggrmap.f90 index adc6dc7b..5a4d44ea 100644 --- a/psb_dgenaggrmap.f90 +++ b/psb_dgenaggrmap.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dgenaggrmap + implicit none integer, intent(in) :: aggr_type type(psb_dspmat_type), intent(in) :: a @@ -56,7 +57,7 @@ subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) if(psb_get_errstatus().ne.0) return info=0 - name = 'psb_bldaggrmat' + name = 'psb_dgenaggrmap' call psb_erractionsave(err_act) ! ! Note. At the time being we are ignoring aggr_type diff --git a/psb_dilu_bld.f90 b/psb_dilu_bld.f90 index 6d1e21c0..c951ea41 100644 --- a/psb_dilu_bld.f90 +++ b/psb_dilu_bld.f90 @@ -49,141 +49,80 @@ !* * !* * !***************************************************************************** -subroutine psb_dilu_bld(a,desc_a,p,upd,info) +subroutine psb_dilu_bld(a,desc_a,p,upd,info,blck) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dilu_bld + implicit none ! ! .. Scalar Arguments .. - integer, intent(out) :: info - ! .. array Arguments .. - type(psb_dspmat_type), intent(in), target :: a - type(psb_dbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a - character, intent(in) :: upd + integer, intent(out) :: info + ! .. array Arguments .. + type(psb_dspmat_type), intent(in), target :: a + type(psb_dbaseprc_type), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + character, intent(in) :: upd + type(psb_dspmat_type), intent(in), optional :: blck ! .. Local Scalars .. integer :: i, j, jj, k, kk, m integer :: int_err(5) character :: trans, unitd - type(psb_dspmat_type) :: blck, atmp real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8 logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false. integer nztota, nztotb, nztmp, nzl, nnr, ir, err_act,& - & n_row, nrow_a,n_col, nhalo, ind, iind, i1,i2,ia + & n_row, nrow_a, ind, iind, i1,i2,ia integer :: ictxt,np,me character(len=20) :: name, ch_err - interface psb_ilu_fct - subroutine psb_dilu_fct(a,l,u,d,info,blck) - use psb_base_mod - 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 - end interface - - interface psb_asmatbld - Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_dspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_dasmatbld - end interface - - interface psb_sp_renum - subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) - use psb_base_mod - use psb_prec_type - implicit none - 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 - if(psb_get_errstatus().ne.0) return info=0 - name='psb_ilu_bld' + name='psb_dilu_bld' call psb_erractionsave(err_act) ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) - m = a%m - if (m < 0) then - info = 10 - int_err(1) = 1 - int_err(2) = m - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif trans = 'N' unitd = 'U' - if (p%iprcparm(n_ovr_) < 0) then - info = 11 - int_err(1) = 1 - int_err(2) = p%iprcparm(n_ovr_) - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif - - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - t1= psb_wtime() - - if(debug) write(0,*)me,': calling psb_asmatbld',p%iprcparm(p_type_),p%iprcparm(n_ovr_) - if (debug) call psb_barrier(ictxt) - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info) - if(info/=0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - t2= psb_wtime() - if (debug) write(0,*)me,': out of psb_asmatbld' - if (debug) call psb_barrier(ictxt) if (allocated(p%av)) then if (size(p%av) < bp_ilu_avsz) then - call psb_errpush(4010,name,a_err='Insufficient av size') - goto 9999 + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) endif - else - call psb_errpush(4010,name,a_err='AV not associated') - goto 9999 + end if + if (.not.allocated(p%av)) then + allocate(p%av(max_avsz),stat=info) + if (info /= 0) then + call psb_errpush(4000,name) + goto 9999 + end if endif !!$ call psb_csprt(50+me,a,head='% (A)') nrow_a = psb_cd_get_local_rows(desc_a) nztota = psb_sp_get_nnzeros(a) - nztotb = psb_sp_get_nnzeros(blck) - if (debug) write(0,*)me,': out get_nnzeros',nztota + if (present(blck)) then + nztota = nztota + psb_sp_get_nnzeros(blck) + end if + if (debug) write(0,*)me,': out get_nnzeros',nztota,a%m,a%k if (debug) call psb_barrier(ictxt) - n_col = psb_cd_get_local_cols(desc_a) - nhalo = n_col-nrow_a n_row = p%desc_data%matrix_data(psb_n_row_) p%av(l_pr_)%m = n_row p%av(l_pr_)%k = n_row p%av(u_pr_)%m = n_row p%av(u_pr_)%k = n_row - call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota+nztotb,info) - if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota+nztotb,info) + call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota,info) + if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota,info) if(info/=0) then info=4010 ch_err='psb_sp_all' @@ -206,100 +145,17 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info) endif - if (debug) then - write(0,*) me,'Done psb_asmatbld' - call psb_barrier(ictxt) - endif - - - if (p%iprcparm(iren_) > 0) then - - ! - ! Here we allocate a full copy to hold local A and received BLK - ! - - nztota = psb_sp_get_nnzeros(a) - nztotb = psb_sp_get_nnzeros(blck) - call psb_sp_all(atmp,nztota+nztotb,info) - if(info/=0) then - info=4011 - call psb_errpush(info,name) - goto 9999 - end if - - - ! write(0,*) 'ILU_BLD ',nztota,nztotb,a%m - - call psb_sp_renum(a,desc_a,blck,p,atmp,info) - - if(info/=0) then - info=4010 - ch_err='psb_sp_renum' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - t3 = psb_wtime() - if (debugprt) then - call psb_barrier(ictxt) - open(40+me) - call psb_csprt(40+me,atmp,head='% Local matrix') - close(40+me) - endif - if (debug) write(0,*) me,' Factoring rows ',& - &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 - - ! - ! Ok, factor the matrix. - ! - t5 = psb_wtime() - blck%m=0 - call psb_ilu_fct(atmp,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_ilu_fct' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_sp_free(atmp,info) - if(info/=0) then - info=4010 - ch_err='psb_sp_free' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - else if (p%iprcparm(iren_) == 0) then - t3 = psb_wtime() - ! This is where we have mo renumbering, thus no need - ! for ATMP - - if (debugprt) then - open(40+me) - call psb_barrier(ictxt) - call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,& - & head='% Local matrix') - if (p%iprcparm(p_type_)==asm_) then - call psb_csprt(40+me,blck,iv=p%desc_data%loc_to_glob,& - & irs=a%m,head='% Received rows') - endif - close(40+me) - endif - - t5= psb_wtime() - if (debug) write(0,*) me,' Going for ilu_fct' - if (debug) call psb_barrier(ictxt) - call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_ilu_fct' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (debug) write(0,*) me,' Done dilu_fct' - endif + ! + ! Ok, factor the matrix. + ! + t5 = psb_wtime() + call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + if(info/=0) then + info=4010 + ch_err='psb_ilu_fct' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if if (debugprt) then @@ -327,28 +183,14 @@ subroutine psb_dilu_bld(a,desc_a,p,upd,info) ! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5 ! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5 - call psb_sp_free(blck,info) - if(info/=0) then - info=4010 - ch_err='psb_sp_free' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (psb_sp_getifld(psb_upd_,p%av(u_pr_),info) /= psb_upd_perm_) then call psb_sp_trimsize(p%av(u_pr_),i1,i2,ia,info) if (info == 0) call psb_sp_reall(p%av(u_pr_),i1,i2,ia,info) - if (info /=0) then - write(0,*) 'Error from trimsize 1',info - endif endif if (psb_sp_getifld(psb_upd_,p%av(l_pr_),info) /= psb_upd_perm_) then call psb_sp_trimsize(p%av(l_pr_),i1,i2,ia,info) if (info == 0) call psb_sp_reall(p%av(l_pr_),i1,i2,ia,info) - if (info /=0) then - write(0,*) 'Error from trimsize 2',info - endif endif diff --git a/psb_dilu_fct.f90 b/psb_dilu_fct.f90 index 950a9b1c..1228d54e 100644 --- a/psb_dilu_fct.f90 +++ b/psb_dilu_fct.f90 @@ -57,7 +57,8 @@ subroutine psb_dilu_fct(a,l,u,d,info,blck) type(psb_dspmat_type), pointer :: blck_ character(len=20) :: name, ch_err logical, parameter :: debug=.false. - name='psb_dcsrlu' + + name='psb_dilu_fct' info = 0 call psb_erractionsave(err_act) ! .. Executable Statements .. diff --git a/psb_dmlprc_aply.f90 b/psb_dmlprc_aply.f90 index 3f675b07..e7652bf8 100644 --- a/psb_dmlprc_aply.f90 +++ b/psb_dmlprc_aply.f90 @@ -84,7 +84,7 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dmlprc_aply implicit none type(psb_desc_type),intent(in) :: desc_data @@ -112,21 +112,7 @@ subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) end type psb_mlprec_wrk_type type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbaseprc_aply - end interface - - name='psb_mlprc_aply' + name='psb_dmlprc_aply' info = 0 call psb_erractionsave(err_act) diff --git a/psb_dmlprc_bld.f90 b/psb_dmlprc_bld.f90 index d26582a8..7037ed38 100644 --- a/psb_dmlprc_bld.f90 +++ b/psb_dmlprc_bld.f90 @@ -37,7 +37,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dmlprc_bld implicit none type(psb_dspmat_type), intent(in), target :: a @@ -51,48 +51,9 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info) character(len=20) :: name, ch_err logical, parameter :: debug=.false. type(psb_dspmat_type) :: ac - - interface psb_baseprc_bld - subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd) - use psb_base_mod - 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 - subroutine psb_dgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) - use psb_base_mod - use psb_prec_type - implicit none - integer, intent(in) :: aggr_type - type(psb_dspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer, allocatable :: ilaggr(:),nlaggr(:) - integer, intent(out) :: info - end subroutine psb_dgenaggrmap - end interface - - interface psb_bldaggrmat - subroutine psb_dbldaggrmat(a,desc_a,ac,desc_ac,p,info) - use psb_base_mod - use psb_prec_type - 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_desc_type), intent(inout) :: desc_ac - type(psb_dbaseprc_type), intent(inout), target :: p - integer, intent(out) :: info - end subroutine psb_dbldaggrmat - end interface - integer :: ictxt, np, me - name='psb_mlprec_bld' + name='psb_dmlprec_bld' if (psb_get_errstatus().ne.0) return info = 0 ictxt = psb_cd_get_context(desc_a) diff --git a/psb_dprc_aply.f90 b/psb_dprc_aply.f90 index 89b3a2c6..76cefd0b 100644 --- a/psb_dprc_aply.f90 +++ b/psb_dprc_aply.f90 @@ -37,7 +37,8 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dprc_aply + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -54,35 +55,7 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name - interface psb_baseprc_aply - subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - real(kind(0.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - real(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_dbaseprc_aply - end interface - - interface psb_mlprc_aply - subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_dbaseprc_type), intent(in) :: baseprecv(:) - real(kind(0.d0)),intent(in) :: alpha,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_dprc_aply' + name='psb_dprec_aply' info = 0 call psb_erractionsave(err_act) @@ -181,7 +154,8 @@ end subroutine psb_dprc_aply subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dprc_aply1 + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -191,27 +165,13 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) character(len=1), optional :: trans logical,parameter :: debug=.false., debugprt=.false. - interface - subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) - use psb_base_mod - use psb_prec_type - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(kind(0.d0)),intent(inout) :: x(:), y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(kind(0.d0)), optional, target :: work(:) - end subroutine psb_dprc_aply - end interface - ! Local variables character :: trans_ integer :: ictxt,np,me,i, err_act real(kind(1.d0)), pointer :: WW(:), w1(:) character(len=20) :: name - name='psb_dprec1' + + name='psb_dprec_aply1' info = 0 call psb_erractionsave(err_act) diff --git a/psb_dprecbld.f90 b/psb_dprecbld.f90 index 7032b6c9..295f5ebc 100644 --- a/psb_dprecbld.f90 +++ b/psb_dprecbld.f90 @@ -1,3 +1,4 @@ + !!$ !!$ !!$ MD2P4 @@ -37,8 +38,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) use psb_base_mod - use psb_prec_type - use psb_prec_mod + use psb_prec_mod, protect => psb_dprecbld Implicit None type(psb_dspmat_type), target :: a @@ -60,7 +60,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_precbld' + name = 'psb_dprecbld' if (debug) write(0,*) 'Entering precbld',P%prec,desc_a%matrix_data(:) info = 0 diff --git a/psb_dprecfree.f90 b/psb_dprecfree.f90 index c1497caf..d8a06b91 100644 --- a/psb_dprecfree.f90 +++ b/psb_dprecfree.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_dprecfree(p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dprecfree + implicit none type(psb_dprec_type), intent(inout) :: p integer, intent(out) :: info @@ -47,7 +48,7 @@ subroutine psb_dprecfree(p,info) if(psb_get_errstatus().ne.0) return info=0 - name = 'psdprecfree' + name = 'psb_dprecfree' call psb_erractionsave(err_act) me=-1 diff --git a/psb_dprecset.f90 b/psb_dprecset.f90 index 5a3ca418..38467236 100644 --- a/psb_dprecset.f90 +++ b/psb_dprecset.f90 @@ -37,7 +37,8 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dprecset + implicit none type(psb_dprec_type), intent(inout) :: p character(len=*), intent(in) :: ptype diff --git a/psb_dslu_bld.f90 b/psb_dslu_bld.f90 index daea1dfa..d579095e 100644 --- a/psb_dslu_bld.f90 +++ b/psb_dslu_bld.f90 @@ -36,142 +36,46 @@ !!$ subroutine psb_dslu_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dslu_bld implicit none - type(psb_dspmat_type), intent(inout) :: a + 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 - character(len=5) :: fmt - character :: upd='F' integer :: i,j,nza,nzb,nzt,ictxt,me,np,err_act logical, parameter :: debug=.false. character(len=20) :: name, ch_err - interface psb_asmatbld - Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_dspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_dasmatbld - end interface - if(psb_get_errstatus().ne.0) return info=0 - name='psb_slu_bld' + name='psb_dslu_bld' call psb_erractionsave(err_act) - ictxt = desc_a%matrix_data(psb_ctxt_) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) - fmt = 'COO' - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - atmp%fida='COO' - if (Debug) then - write(0,*) me, 'SPLUBLD: Calling csdp' - call psb_barrier(ictxt) - endif - - call psb_csdp(a,atmp,info) - if(info /= 0) then - info=4010 - ch_err='psb_csdp' - call psb_errpush(info,name,a_err=ch_err) + if (toupper(a%fida) /= 'CSR') then + write(0,*) 'Unimplemented input to SLU_BLD' goto 9999 - end if - nza = atmp%infoa(psb_nnz_) - if (Debug) then - write(0,*) me, 'SPLUBLD: Done csdp',info,nza,atmp%m,atmp%k - call psb_barrier(ictxt) endif - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=fmt) - if(info /= 0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - nzb = blck%infoa(psb_nnz_) - if (Debug) then - write(0,*) me, 'SPLUBLD: Done asmatbld',info,nzb,blck%fida - call psb_barrier(ictxt) - endif - if (nzb > 0 ) then - if (size(atmp%aspk) psb_dsp_renum + implicit none ! .. array Arguments .. @@ -57,24 +58,32 @@ subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) if (psb_get_errstatus().ne.0) return info=0 - name='apply_renum' + name='psb_dsp_renum' call psb_erractionsave(err_act) ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) -!!!!!!!!!!!!!!!! CHANGE FOR NON-CSR A + ! + ! CHANGE: Start with a COO atmp. Then change if/when necessary. + ! Exit with a COO atmp. ! ! Renumbering type: ! 1. Global column indices ! (2. GPS band reduction disabled for the time being) + nztota=psb_sp_get_nnzeros(a) + nztotb=psb_sp_get_nnzeros(blck) + call psb_sp_reall(atmp,nztota+nztotb,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_sp_reall') + goto 9999 + end if + atmp%fida='COO' + call psb_csdp(a,atmp,info) + call psb_rwextd(a%m+blck%m,atmp,info,blck,rowscale=.false.) if (p%iprcparm(iren_)==renum_glb_) then - atmp%m = a%m + blck%m - atmp%k = a%k - atmp%fida='CSR' - atmp%descra = 'GUN' ! This is the renumbering coherent with global indices.. mglob = psb_cd_get_global_rows(desc_a) @@ -104,112 +113,10 @@ subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) enddo t3 = psb_wtime() - ! Build ATMP with new numbering. - nztmp=size(atmp%aspk) - allocate(itmp(max(8,atmp%m+2,nztmp+2)),rtmp(atmp%m),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - j = 1 - atmp%ia2(1) = 1 - do i=1, atmp%m - ir = p%perm(i) - - if (ir <= a%m ) then - - nzl = a%ia2(ir+1) - a%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = a%ia2(ir) - k=0 - do kk=1, nzl - if (a%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = a%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) - endif - enddo - call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else if (ir <= atmp%m ) then - - ir = ir - a%m - nzl = blck%ia2(ir+1) - blck%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = blck%ia2(ir) - k=0 - do kk=1, nzl - if (blck%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = blck%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) - endif - enddo - call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else - write(0,*) 'Row index error 1 :',i,ir - endif - - j = j + k - atmp%ia2(i+1) = j - - enddo - - t4 = psb_wtime() - - - deallocate(itmp,itmp2,rtmp) - else if (p%iprcparm(iren_)==renum_gps_) then - - atmp%m = a%m + blck%m - atmp%k = a%k - atmp%fida='CSR' - atmp%descra = 'GUN' - do i=1, a%m - atmp%ia2(i) = a%ia2(i) - do j= a%ia2(i), a%ia2(i+1)-1 - atmp%ia1(j) = a%ia1(j) - enddo - enddo - atmp%ia2(a%m+1) = a%ia2(a%m+1) - nztota = atmp%ia2(a%m+1) -1 - if (blck%m>0) then - do i=1, blck%m - atmp%ia2(a%m+i) = nztota+blck%ia2(i) - do j= blck%ia2(i), blck%ia2(i+1)-1 - atmp%ia1(nztota+j) = blck%ia1(j) - enddo - enddo - atmp%ia2(atmp%m+1) = nztota+blck%ia2(blck%m+1) - endif - nztmp = atmp%ia2(atmp%m+1) - 1 - - + + call psb_ipcoo2csr(atmp,info) + nztmp = psb_sp_get_nnzeros(atmp) ! This is a renumbering with Gibbs-Poole-Stockmeyer ! band reduction. Switched off for now. To be fixed, ! gps_reduction should get p%perm. @@ -263,90 +170,26 @@ subroutine psb_dsp_renum(a,desc_a,blck,p,atmp,info) p%invperm(p%perm(k)) = k enddo t3 = psb_wtime() - - ! Build ATMP with new numbering. - - allocate(itmp2(max(8,atmp%m+2,nztmp+2)),rtmp(atmp%m),stat=info) - if (info /= 0) then - call psb_errpush(4010,name,a_err='Allocate') - goto 9999 - end if - - j = 1 - atmp%ia2(1) = 1 - do i=1, atmp%m - ir = p%perm(i) - - if (ir <= a%m ) then - - nzl = a%ia2(ir+1) - a%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = a%ia2(ir) - k=0 - do kk=1, nzl - if (a%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = a%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(a%ia1(jj+kk-1)) - endif - enddo - call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else if (ir <= atmp%m ) then - - ir = ir - a%m - nzl = blck%ia2(ir+1) - blck%ia2(ir) - if (nzl > size(rtmp)) then - call psb_realloc(nzl,rtmp,info) - if(info/=0) then - info=4010 - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - endif - jj = blck%ia2(ir) - k=0 - do kk=1, nzl - if (blck%ia1(jj+kk-1)<=atmp%m) then - k = k + 1 - rtmp(k) = blck%aspk(jj+kk-1) - atmp%ia1(j+k-1) = p%invperm(blck%ia1(jj+kk-1)) - endif - enddo - call psb_msort(atmp%ia1(j:j+k-1),ix=itmp2) - do kk=1,k - atmp%aspk(j+kk-1) = rtmp(itmp2(kk)) - enddo - - else - write(0,*) 'Row index error 1 :',i,ir - endif - - j = j + k - atmp%ia2(i+1) = j - - enddo - - t4 = psb_wtime() - - - - deallocate(itmp,itmp2,rtmp) + + call psb_ipcsr2coo(atmp,info) end if + ! Rebuild ATMP with new numbering. + + nztmp=psb_sp_get_nnzeros(atmp) + do i=1,nztmp + atmp%ia1(i) = p%perm(a%ia1(i)) + atmp%ia2(i) = p%invperm(a%ia2(i)) + end do + call psb_fixcoo(atmp,info) + if (info /= 0) then + call psb_errpush(4010,name,a_err='psb_fixcoo') + goto 9999 + end if + + t4 = psb_wtime() + call psb_erractionrestore(err_act) return diff --git a/psb_dumf_bld.f90 b/psb_dumf_bld.f90 index 1569f979..88573d90 100644 --- a/psb_dumf_bld.f90 +++ b/psb_dumf_bld.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_dumf_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_dumf_bld + implicit none type(psb_dspmat_type), intent(inout) :: a @@ -45,140 +46,43 @@ subroutine psb_dumf_bld(a,desc_a,p,info) integer, intent(out) :: info - type(psb_dspmat_type) :: blck, atmp - character(len=5) :: fmt - character :: upd='F' integer :: i,j,nza,nzb,nzt,ictxt,me,np,err_act integer :: i_err(5) logical, parameter :: debug=.false. character(len=20) :: name, ch_err - interface psb_asmatbld - Subroutine psb_dasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_dspmat_type), Intent(in) :: a - Type(psb_dspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_dasmatbld - end interface - info=0 - name='psb_umf_bld' + name='psb_dumf_bld' call psb_erractionsave(err_act) ictxt = desc_A%matrix_data(psb_ctxt_) call psb_info(ictxt, me, np) - fmt = 'COO' - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - atmp%fida='COO' - if (Debug) then - write(0,*) me, 'UMFBLD: Calling csdp' - call psb_barrier(ictxt) - endif - - call psb_dcsdp(a,atmp,info) - if(info /= 0) then - info=4010 - ch_err='psb_dcsdp' - call psb_errpush(info,name,a_err=ch_err) + if (toupper(a%fida) /= 'CSC') then + write(0,*) 'Unimplemented input to UMF_BLD' goto 9999 - end if - nza = psb_sp_get_nnzeros(atmp) - nzb = psb_sp_get_nnzeros(a) - - if (Debug) then - write(0,*) me, 'UMFBLD: Done csdp',info,nza,atmp%m,atmp%k,nzb - call psb_barrier(ictxt) endif - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=fmt) - if(info /= 0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - nzb = psb_sp_get_nnzeros(blck) - if (Debug) then - write(0,*) me, 'UMFBLD: Done asmatbld',info,nzb,blck%fida - call psb_barrier(ictxt) - endif - if (nzb > 0 ) then - if (size(atmp%aspk)=psb_none_).and.(ip<=psb_square_root_)) return end function is_legal_prolong function is_legal_restrict(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_restrict is_legal_restrict = ((ip==psb_nohalo_).or.(ip==psb_halo_)) return end function is_legal_restrict function is_legal_ml_type(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_type @@ -565,6 +577,7 @@ contains return end function is_legal_ml_type function is_legal_ml_aggr_kind(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_aggr_kind @@ -572,6 +585,7 @@ contains return end function is_legal_ml_aggr_kind function is_legal_ml_smooth_pos(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_smooth_pos @@ -579,6 +593,7 @@ contains return end function is_legal_ml_smooth_pos function is_legal_ml_smth_kind(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_smth_kind @@ -586,6 +601,7 @@ contains return end function is_legal_ml_smth_kind function is_legal_ml_coarse_mat(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_coarse_mat @@ -593,6 +609,7 @@ contains return end function is_legal_ml_coarse_mat function is_legal_ml_fact(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_fact @@ -600,6 +617,7 @@ contains return end function is_legal_ml_fact function is_legal_ml_lev(ip) + use psb_base_mod integer, intent(in) :: ip logical :: is_legal_ml_lev @@ -607,6 +625,7 @@ contains return end function is_legal_ml_lev function is_legal_omega(ip) + use psb_base_mod real(kind(1.d0)), intent(in) :: ip logical :: is_legal_omega @@ -614,6 +633,7 @@ contains return end function is_legal_omega function is_legal_ml_eps(ip) + use psb_base_mod real(kind(1.d0)), intent(in) :: ip logical :: is_legal_ml_eps @@ -623,6 +643,7 @@ contains subroutine psb_icheck_def(ip,name,id,is_legal) + use psb_base_mod integer, intent(inout) :: ip integer, intent(in) :: id character(len=*), intent(in) :: name @@ -640,6 +661,7 @@ contains end subroutine psb_icheck_def subroutine psb_dcheck_def(ip,name,id,is_legal) + use psb_base_mod real(kind(1.d0)), intent(inout) :: ip real(kind(1.d0)), intent(in) :: id character(len=*), intent(in) :: name @@ -823,6 +845,7 @@ contains function pr_to_str(iprec) + use psb_base_mod integer, intent(in) :: iprec character(len=10) :: pr_to_str diff --git a/psb_zasmatbld.f90 b/psb_zasmatbld.f90 index 822ee0b8..f112be2c 100644 --- a/psb_zasmatbld.f90 +++ b/psb_zasmatbld.f90 @@ -54,7 +54,8 @@ Subroutine psb_zasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zasmatbld + Implicit None ! .. Array Arguments .. @@ -77,6 +78,7 @@ Subroutine psb_zasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) & tot_recv, ircode, n_row,nhalo, nrow_a,err_act Logical,Parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err + name='psb_zasmatbld' if(psb_get_errstatus().ne.0) return info=0 diff --git a/psb_zbaseprc_aply.f90 b/psb_zbaseprc_aply.f90 index d683df05..2032f797 100644 --- a/psb_zbaseprc_aply.f90 +++ b/psb_zbaseprc_aply.f90 @@ -40,7 +40,8 @@ subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! where K is a a basic preconditioner stored in prec ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zbaseprc_aply + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -59,20 +60,6 @@ subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7 logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err - - interface psb_bjac_aply - subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type), intent(in) :: desc_data - type(psb_zbaseprc_type), intent(in) :: prec - complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - complex(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_zbjac_aply - end interface name='psb_zbaseprc_aply' info = 0 diff --git a/psb_zbaseprc_bld.f90 b/psb_zbaseprc_bld.f90 index 1b0abc28..d5e51df1 100644 --- a/psb_zbaseprc_bld.f90 +++ b/psb_zbaseprc_bld.f90 @@ -37,7 +37,8 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zbaseprc_bld + Implicit None type(psb_zspmat_type), target :: a @@ -46,52 +47,6 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) integer, intent(out) :: info character, intent(in), optional :: upd - interface psb_diagsc_bld - subroutine psb_zdiagsc_bld(a,desc_data,p,upd,info) - use psb_base_mod - use psb_prec_type - integer, intent(out) :: info - type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_zdiagsc_bld - end interface - - interface psb_ilu_bld - subroutine psb_zilu_bld(a,desc_data,p,upd,info) - use psb_base_mod - use psb_prec_type - integer, intent(out) :: info - type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(inout) :: p - character, intent(in) :: upd - end subroutine psb_zilu_bld - end interface - - interface psb_slu_bld - subroutine psb_zslu_bld(a,desc_a,p,info) - use psb_base_mod - use psb_prec_type - type(psb_zspmat_type), intent(inout) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_zbaseprc_type), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_zslu_bld - end interface - - interface psb_umf_bld - subroutine psb_zumf_bld(a,desc_a,p,info) - use psb_base_mod - use psb_prec_type - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_zbaseprc_type), intent(inout) :: p - integer, intent(out) :: info - end subroutine psb_zumf_bld - end interface - ! Local scalars Integer :: err, nnzero, n_row, n_col,I,j,k,ictxt,& & me,mycol,np,npcol,mglob,lw, mtype, nrg, nzg, err_act @@ -108,7 +63,7 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_baseprc_bld' + name = 'psb_zbaseprc_bld' if (debug) write(0,*) 'Entering baseprc_bld' info = 0 @@ -140,13 +95,8 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) call psb_check_def(p%iprcparm(p_type_),'base_prec',& & diag_,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) + + call psb_nullify_desc(p%desc_data) select case(p%iprcparm(p_type_)) case (noprec_) @@ -183,59 +133,16 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) 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' + if (debug) write(0,*)me, ': Calling PSB_BJAC_BLD' if (debug) call psb_barrier(ictxt) - 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 (debug) call psb_barrier(ictxt) - 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(debug) write(0,*)me,': Done umf_bld ',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??' + call psb_bjac_bld(a,desc_a,p,iupd,info) + if(info /= 0) then info=4010 - ch_err='Inconsistent prec f_none_' - call psb_errpush(info,name,a_err=ch_err) + call psb_errpush(info,name,a_err='psb_bjac_bld') goto 9999 + end if - case default - write(0,*) 'Unknown factor type in baseprc_bld bja/asm: ',& - &p%iprcparm(f_type_) - info=4010 - ch_err='Unknown f_type_' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end select case default info=4010 ch_err='Unknown p_type_' diff --git a/psb_zbjac_aply.f90 b/psb_zbjac_aply.f90 index f92d589a..422bca17 100644 --- a/psb_zbjac_aply.f90 +++ b/psb_zbjac_aply.f90 @@ -43,7 +43,8 @@ subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zbjac_aply + implicit none type(psb_desc_type), intent(in) :: desc_data @@ -63,7 +64,7 @@ subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name, ch_err - name='psb_bjac_aply' + name='psb_zbjac_aply' info = 0 call psb_erractionsave(err_act) diff --git a/psb_zbjac_bld.f90 b/psb_zbjac_bld.f90 new file mode 100644 index 00000000..fa572dd0 --- /dev/null +++ b/psb_zbjac_bld.f90 @@ -0,0 +1,374 @@ +!!$ +!!$ +!!$ MD2P4 +!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS +!!$ 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 Second 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 MD2P4 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 MD2P4 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. +!!$ +!!$ +!***************************************************************************** +!* * +!* This is where the action takes place. * +!* ASMATBLD does the setup: building the prec descriptor plus retrieving * +!* matrix rows if needed * +!* * +!* * +!* * +!* * +!* some open code does the renumbering * +!* * +!* * +!* * +!* * +!***************************************************************************** +subroutine psb_zbjac_bld(a,desc_a,p,upd,info) + use psb_base_mod + use psb_prec_mod, mld_protect_name => psb_zbjac_bld + + implicit none + ! + ! .. Scalar Arguments .. + integer, intent(out) :: info + ! .. array Arguments .. + type(psb_zspmat_type), intent(in), target :: a + type(psb_zbaseprc_type), intent(inout) :: p + type(psb_desc_type), intent(in) :: desc_a + character, intent(in) :: upd + + ! .. Local Scalars .. + integer :: i, j, jj, k, kk, m + integer :: int_err(5) + character :: trans, unitd + type(psb_zspmat_type) :: blck, atmp + real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8 + logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false. + integer nztota, nztotb, nztmp, nzl, nnr, ir, err_act,& + & n_row, nrow_a,n_col, nhalo, ind, iind, i1,i2,ia + integer :: ictxt,np,me + character(len=20) :: name, ch_err + character(len=5), parameter :: coofmt='COO' + + if(psb_get_errstatus().ne.0) return + info=0 + name='psb_bjac_bld' + call psb_erractionsave(err_act) + + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt, me, np) + + m = a%m + if (m < 0) then + info = 10 + int_err(1) = 1 + int_err(2) = m + call psb_errpush(info,name,i_err=int_err) + goto 9999 + endif + trans = 'N' + unitd = 'U' + if (p%iprcparm(n_ovr_) < 0) then + info = 11 + int_err(1) = 1 + int_err(2) = p%iprcparm(n_ovr_) + call psb_errpush(info,name,i_err=int_err) + goto 9999 + endif + + call psb_nullify_sp(blck) + call psb_nullify_sp(atmp) + + t1= psb_wtime() + + if(debug) write(0,*)me,': calling psb_asmatbld',p%iprcparm(p_type_),p%iprcparm(n_ovr_) + if (debug) call psb_barrier(ictxt) + call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& + & blck,desc_a,upd,p%desc_data,info,outfmt=coofmt) + + if(info/=0) then + info=4010 + ch_err='psb_asmatbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + t2= psb_wtime() + if (debug) write(0,*)me,': out of psb_asmatbld' + if (debug) call psb_barrier(ictxt) + + + if (p%iprcparm(iren_) > 0) then + + ! + ! Here we allocate a full copy to hold local A and received BLK + ! Done inside sp_renum. + ! + + call psb_sp_renum(a,desc_a,blck,p,atmp,info) + + if(info/=0) then + info=4010 + ch_err='psb_sp_renum' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + t3 = psb_wtime() + if (debugprt) then + call psb_barrier(ictxt) + open(40+me) + call psb_csprt(40+me,atmp,head='% Local matrix') + close(40+me) + endif + if (debug) write(0,*) me,' Factoring rows ',& + &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 + + select case(p%iprcparm(f_type_)) + + case(f_ilu_n_,f_ilu_e_) + + call psb_ipcoo2csr(atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + call psb_ilu_bld(atmp,p%desc_data,p,upd,info) + + if(info/=0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ilu_bld') + goto 9999 + end if + + if (debugprt) then + + open(80+me) + + call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') + write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m + do i=1,p%av(l_pr_)%m + write(80+me,*) i,i,p%d(i) + enddo + call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') + + close(80+me) + endif + + + case(f_slu_) + + call psb_ipcoo2csr(atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ipcoo2csr') + goto 9999 + end if + + call psb_slu_bld(atmp,p%desc_data,p,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='slu_bld') + goto 9999 + end if + + case(f_umf_) + + call psb_ipcoo2csc(atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ipcoo2csc') + goto 9999 + end if + + call psb_umf_bld(atmp,p%desc_data,p,info) + if(debug) write(0,*)me,': Done umf_bld ',info + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='umf_bld') + goto 9999 + end if + + case(f_none_) + info=4010 + call psb_errpush(info,name,a_err='Inconsistent prec f_none_') + goto 9999 + + case default + info=4010 + call psb_errpush(info,name,a_err='Unknown f_type_') + goto 9999 + end select + + + + call psb_sp_free(atmp,info) + if(info/=0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_free') + goto 9999 + end if + + + else if (p%iprcparm(iren_) == 0) then + + + select case(p%iprcparm(f_type_)) + + case(f_ilu_n_,f_ilu_e_) + + call psb_ipcoo2csr(blck,info,rwshr=.true.) + + if (info==0) call psb_ilu_bld(a,desc_a,p,upd,info,blck=blck) + + if(info/=0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ilu_bld') + goto 9999 + end if + + if (debugprt) then + + open(80+me) + + call psb_csprt(80+me,p%av(l_pr_),head='% Local L factor') + write(80+me,*) '% Diagonal: ',p%av(l_pr_)%m + do i=1,p%av(l_pr_)%m + write(80+me,*) i,i,p%d(i) + enddo + call psb_csprt(80+me,p%av(u_pr_),head='% Local U factor') + + close(80+me) + endif + + + + case(f_slu_) + + atmp%fida='COO' + call psb_csdp(a,atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_csdp') + goto 9999 + end if + call psb_rwextd(atmp%m+blck%m,atmp,info,blck,rowscale=.true.) + call psb_ipcoo2csr(atmp,info) + call psb_slu_bld(atmp,p%desc_data,p,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='slu_bld') + goto 9999 + end if + + call psb_sp_free(atmp,info) + if(info/=0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_free') + goto 9999 + end if + + + case(f_umf_) + + atmp%fida='COO' + atmp%fida='COO' + call psb_csdp(a,atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_csdp') + goto 9999 + end if + call psb_rwextd(atmp%m+blck%m,atmp,info,blck,rowscale=.true.) + if (info == 0) call psb_ipcoo2csc(atmp,info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_ipcoo2csc') + goto 9999 + end if + + call psb_umf_bld(atmp,p%desc_data,p,info) + if(debug) write(0,*)me,': Done umf_bld ',info + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err='umf_bld') + goto 9999 + end if + + call psb_sp_free(atmp,info) + if(info/=0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_free') + goto 9999 + end if + + + case(f_none_) + info=4010 + call psb_errpush(info,name,a_err='Inconsistent prec f_none_') + goto 9999 + + case default + info=4010 + call psb_errpush(info,name,a_err='Unknown f_type_') + goto 9999 + end select + + + endif + + t6 = psb_wtime() + + call psb_sp_free(blck,info) + if(info/=0) then + info=4010 + ch_err='psb_sp_free' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if (debug) write(0,*) me,'End of ilu_bld' + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + + +end subroutine psb_zbjac_bld + + diff --git a/psb_zbldaggrmat.f90 b/psb_zbldaggrmat.f90 index 3d9d2c5d..0e9347b8 100644 --- a/psb_zbldaggrmat.f90 +++ b/psb_zbldaggrmat.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zbldaggrmat + implicit none type(psb_zspmat_type), intent(in), target :: a @@ -49,6 +50,7 @@ subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info) logical, parameter :: aggr_dump=.false. integer ::ictxt,np,me, err_act character(len=20) :: name, ch_err + name='psb_zbldaggrmat' if(psb_get_errstatus().ne.0) return info=0 diff --git a/psb_zdiagsc_bld.f90 b/psb_zdiagsc_bld.f90 index 077d3179..bd37d9e0 100644 --- a/psb_zdiagsc_bld.f90 +++ b/psb_zdiagsc_bld.f90 @@ -37,7 +37,8 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zdiagsc_bld + Implicit None type(psb_zspmat_type), target :: a @@ -62,7 +63,7 @@ subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_diagsc_bld' + name = 'psb_zdiagsc_bld' if (debug) write(0,*) 'Entering diagsc_bld' info = 0 diff --git a/psb_zgenaggrmap.f90 b/psb_zgenaggrmap.f90 index c203616a..00b81c07 100644 --- a/psb_zgenaggrmap.f90 +++ b/psb_zgenaggrmap.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_zgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zgenaggrmap + implicit none integer, intent(in) :: aggr_type type(psb_zspmat_type), intent(in) :: a @@ -56,7 +57,7 @@ subroutine psb_zgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) if(psb_get_errstatus().ne.0) return info=0 - name = 'psb_bldaggrmat' + name = 'psb_zgenaggrmap' call psb_erractionsave(err_act) ! ! Note. At the time being we are ignoring aggr_type diff --git a/psb_zilu_bld.f90 b/psb_zilu_bld.f90 index 6a14bb2d..44f31d51 100644 --- a/psb_zilu_bld.f90 +++ b/psb_zilu_bld.f90 @@ -49,9 +49,10 @@ !* * !* * !***************************************************************************** -subroutine psb_zilu_bld(a,desc_a,p,upd,info) +subroutine psb_zilu_bld(a,desc_a,p,upd,info,blck) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zilu_bld + implicit none ! ! .. Scalar Arguments .. @@ -61,128 +62,66 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info) type(psb_zbaseprc_type), intent(inout) :: p type(psb_desc_type), intent(in) :: desc_a character, intent(in) :: upd - + type(psb_zspmat_type), intent(in), optional :: blck ! .. Local Scalars .. integer :: i, j, jj, k, kk, m integer :: int_err(5) character :: trans, unitd - type(psb_zspmat_type) :: blck, atmp real(kind(1.d0)) :: t1,t2,t3,t4,t5,t6, t7, t8 logical, parameter :: debugprt=.false., debug=.false., aggr_dump=.false. integer nztota, nztotb, nztmp, nzl, nnr, ir, err_act,& - & n_row, nrow_a,n_col, nhalo, ind, iind, i1,i2,ia + & n_row, nrow_a, ind, iind, i1,i2,ia integer :: ictxt,np,me character(len=20) :: name, ch_err - interface psb_ilu_fct - subroutine psb_zilu_fct(a,l,u,d,info,blck) - use psb_base_mod - use psb_prec_type - integer, intent(out) :: info - type(psb_zspmat_type),intent(in) :: a - type(psb_zspmat_type),intent(inout) :: l,u - type(psb_zspmat_type),intent(in), optional, target :: blck - complex(kind(1.d0)), intent(inout) :: d(:) - end subroutine psb_zilu_fct - end interface - - interface psb_asmatbld - Subroutine psb_zasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_zspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_zasmatbld - end interface - - interface psb_sp_renum - subroutine psb_zsp_renum(a,desc_a,blck,p,atmp,info) - use psb_base_mod - use psb_prec_type - type(psb_zspmat_type), intent(in) :: a,blck - type(psb_zspmat_type), intent(inout) :: atmp - type(psb_zbaseprc_type), intent(inout) :: p - type(psb_desc_type), intent(in) :: desc_a - integer, intent(out) :: info - end subroutine psb_zsp_renum - end interface - if(psb_get_errstatus().ne.0) return info=0 - name='psb_ilu_bld' + name='psb_zilu_bld' call psb_erractionsave(err_act) ictxt=psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) - m = a%m - if (m < 0) then - info = 10 - int_err(1) = 1 - int_err(2) = m - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif trans = 'N' unitd = 'U' - if (p%iprcparm(n_ovr_) < 0) then - info = 11 - int_err(1) = 1 - int_err(2) = p%iprcparm(n_ovr_) - call psb_errpush(info,name,i_err=int_err) - goto 9999 - endif - - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - t1= psb_wtime() - - if(debug) write(0,*)me,': calling psb_asmatbld',p%iprcparm(p_type_),p%iprcparm(n_ovr_) - if (debug) call psb_barrier(ictxt) - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info) - if(info/=0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - t2= psb_wtime() - if (debug) write(0,*)me,': out of psb_asmatbld' - if (debug) call psb_barrier(ictxt) if (allocated(p%av)) then if (size(p%av) < bp_ilu_avsz) then - call psb_errpush(4010,name,a_err='Insufficient av size') - goto 9999 + do i=1,size(p%av) + call psb_sp_free(p%av(i),info) + if (info /= 0) then + ! Actually, we don't care here about this. + ! Just let it go. + ! return + end if + enddo + deallocate(p%av,stat=info) endif - else - call psb_errpush(4010,name,a_err='AV not associated') - goto 9999 + end if + if (.not.allocated(p%av)) then + allocate(p%av(max_avsz),stat=info) + if (info /= 0) then + call psb_errpush(4000,name) + goto 9999 + end if endif +!!$ call psb_csprt(50+me,a,head='% (A)') nrow_a = psb_cd_get_local_rows(desc_a) nztota = psb_sp_get_nnzeros(a) - nztotb = psb_sp_get_nnzeros(blck) - if (debug) write(0,*)me,': out get_nnzeros',nztota + if (present(blck)) then + nztota = nztota + psb_sp_get_nnzeros(blck) + end if + if (debug) write(0,*)me,': out get_nnzeros',nztota,a%m,a%k if (debug) call psb_barrier(ictxt) - n_col = psb_cd_get_local_cols(desc_a) - nhalo = n_col-nrow_a n_row = p%desc_data%matrix_data(psb_n_row_) p%av(l_pr_)%m = n_row p%av(l_pr_)%k = n_row p%av(u_pr_)%m = n_row p%av(u_pr_)%k = n_row - call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota+nztotb,info) - if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota+nztotb,info) + call psb_sp_all(n_row,n_row,p%av(l_pr_),nztota,info) + if (info == 0) call psb_sp_all(n_row,n_row,p%av(u_pr_),nztota,info) if(info/=0) then info=4010 ch_err='psb_sp_all' @@ -205,101 +144,17 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info) endif - if (debug) then - write(0,*) me,'Done psb_asmatbld' - call psb_barrier(ictxt) - endif - - - if (p%iprcparm(iren_) > 0) then - - ! - ! Here we allocate a full copy to hold local A and received BLK - ! - - nztota = psb_sp_get_nnzeros(a) - nztotb = psb_sp_get_nnzeros(blck) - - call psb_sp_all(atmp,nztota+nztotb,info) - if(info/=0) then - info=4011 - call psb_errpush(info,name) - goto 9999 - end if - - - ! write(0,*) 'ILU_BLD ',nztota,nztotb,a%m - - call psb_sp_renum(a,desc_a,blck,p,atmp,info) - - if(info/=0) then - info=4010 - ch_err='psb_sp_renum' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - t3 = psb_wtime() - if (debugprt) then - call psb_barrier(ictxt) - open(40+me) - call psb_csprt(40+me,atmp,head='% Local matrix') - close(40+me) - endif - if (debug) write(0,*) me,' Factoring rows ',& - &atmp%m,a%m,blck%m,atmp%ia2(atmp%m+1)-1 - - ! - ! Ok, factor the matrix. - ! - t5 = psb_wtime() - blck%m=0 - call psb_ilu_fct(atmp,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_ilu_fct' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - call psb_sp_free(atmp,info) - if(info/=0) then - info=4010 - ch_err='psb_sp_free' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - - else if (p%iprcparm(iren_) == 0) then - t3 = psb_wtime() - ! This is where we have mo renumbering, thus no need - ! for ATMP - - if (debugprt) then - open(40+me) - call psb_barrier(ictxt) - call psb_csprt(40+me,a,iv=p%desc_data%loc_to_glob,& - & head='% Local matrix') - if (p%iprcparm(p_type_)==asm_) then - call psb_csprt(40+me,blck,iv=p%desc_data%loc_to_glob,& - & irs=a%m,head='% Received rows') - endif - close(40+me) - endif - - t5= psb_wtime() - if (debug) write(0,*) me,' Going for ilu_fct' - if (debug) call psb_barrier(ictxt) - call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) - if(info/=0) then - info=4010 - ch_err='psb_ilu_fct' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (debug) write(0,*) me,' Done dilu_fct' - endif + ! + ! Ok, factor the matrix. + ! + t5 = psb_wtime() + call psb_ilu_fct(a,p%av(l_pr_),p%av(u_pr_),p%d,info,blck=blck) + if(info/=0) then + info=4010 + ch_err='psb_ilu_fct' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if if (debugprt) then @@ -318,6 +173,8 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info) close(80+me) endif +!!$ call psb_csprt(60+me,a,head='% (A)') + ! ierr = MPE_Log_event( ifcte, 0, "st SIMPLE" ) t6 = psb_wtime() @@ -325,14 +182,6 @@ subroutine psb_zilu_bld(a,desc_a,p,upd,info) ! write(0,'(i3,1x,a,3(1x,g18.9))') me,'renum/factor time',t3-t2,t6-t5 ! if (me==0) write(0,'(a,3(1x,g18.9))') 'renum/factor time',t3-t2,t6-t5 - call psb_sp_free(blck,info) - if(info/=0) then - info=4010 - ch_err='psb_sp_free' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (psb_sp_getifld(psb_upd_,p%av(u_pr_),info) /= psb_upd_perm_) then call psb_sp_trimsize(p%av(u_pr_),i1,i2,ia,info) if (info == 0) call psb_sp_reall(p%av(u_pr_),i1,i2,ia,info) diff --git a/psb_zilu_fct.f90 b/psb_zilu_fct.f90 index 0c00375f..a3a7f236 100644 --- a/psb_zilu_fct.f90 +++ b/psb_zilu_fct.f90 @@ -42,7 +42,6 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck) ! ! use psb_base_mod - use psb_prec_type implicit none ! .. Scalar Arguments .. integer, intent(out) :: info @@ -55,7 +54,8 @@ subroutine psb_zilu_fct(a,l,u,d,info,blck) integer :: i, j, jj, k, kk, l1, l2, ll, low1, low2,m,ma,err_act type(psb_zspmat_type), pointer :: blck_ character(len=20) :: name, ch_err - name='psb_zcsrlu' + + name='psb_zilu_fct' info = 0 call psb_erractionsave(err_act) ! .. Executable Statements .. diff --git a/psb_zmlprc_aply.f90 b/psb_zmlprc_aply.f90 index 2d09808b..c2d28886 100644 --- a/psb_zmlprc_aply.f90 +++ b/psb_zmlprc_aply.f90 @@ -84,7 +84,8 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) ! use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zmlprc_aply + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -111,21 +112,7 @@ subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) end type psb_mlprec_wrk_type type(psb_mlprec_wrk_type), allocatable :: mlprec_wrk(:) - interface psb_baseprc_aply - subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(in) :: prec - complex(kind(1.d0)),intent(inout) :: x(:), y(:) - complex(kind(1.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - complex(kind(1.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_zbaseprc_aply - end interface - - name='psb_mlprc_aply' + name='psb_zmlprc_aply' info = 0 call psb_erractionsave(err_act) diff --git a/psb_zmlprc_bld.f90 b/psb_zmlprc_bld.f90 index b53bc401..98e6d888 100644 --- a/psb_zmlprc_bld.f90 +++ b/psb_zmlprc_bld.f90 @@ -37,7 +37,8 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zmlprc_bld + implicit none type(psb_zspmat_type), intent(in), target :: a @@ -51,48 +52,9 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info) character(len=20) :: name, ch_err logical, parameter :: debug=.false. type(psb_zspmat_type) :: ac - - interface psb_baseprc_bld - subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd) - use psb_base_mod - use psb_prec_type - type(psb_zspmat_type), target :: a - type(psb_desc_type), intent(in) :: desc_a - type(psb_zbaseprc_type),intent(inout) :: p - integer, intent(out) :: info - character, intent(in), optional :: upd - end subroutine psb_zbaseprc_bld - end interface - - interface psb_genaggrmap - subroutine psb_zgenaggrmap(aggr_type,a,desc_a,nlaggr,ilaggr,info) - use psb_base_mod - use psb_prec_type - implicit none - integer, intent(in) :: aggr_type - type(psb_zspmat_type), intent(in) :: a - type(psb_desc_type), intent(in) :: desc_a - integer, allocatable :: ilaggr(:),nlaggr(:) - integer, intent(out) :: info - end subroutine psb_zgenaggrmap - end interface - - interface psb_bldaggrmat - subroutine psb_zbldaggrmat(a,desc_a,ac,desc_ac,p,info) - use psb_base_mod - use psb_prec_type - type(psb_zspmat_type), intent(in), target :: a - type(psb_zbaseprc_type), intent(inout),target :: p - type(psb_zspmat_type), intent(out),target :: ac - type(psb_desc_type), intent(in) :: desc_a - type(psb_desc_type), intent(inout) :: desc_ac - integer, intent(out) :: info - end subroutine psb_zbldaggrmat - end interface - integer :: ictxt, np, me - name='psb_mlprec_bld' + name='psb_zmlprec_bld' if (psb_get_errstatus().ne.0) return info = 0 ictxt = psb_cd_get_context(desc_a) diff --git a/psb_zprc_aply.f90 b/psb_zprc_aply.f90 index 70ebc1ec..c12d2c96 100644 --- a/psb_zprc_aply.f90 +++ b/psb_zprc_aply.f90 @@ -37,7 +37,8 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zprc_aply + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -53,34 +54,6 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) integer :: ictxt,np,me,err_act logical,parameter :: debug=.false., debugprt=.false. character(len=20) :: name - - interface psb_baseprc_aply - subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(in) :: prec - complex(kind(0.d0)),intent(inout) :: x(:), y(:) - complex(kind(0.d0)),intent(in) :: alpha,beta - character(len=1) :: trans - complex(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_zbaseprc_aply - end interface - - interface psb_mlprc_aply - subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_zbaseprc_type), intent(in) :: baseprecv(:) - complex(kind(0.d0)),intent(in) :: alpha,beta - complex(kind(0.d0)),intent(inout) :: x(:), y(:) - character :: trans - complex(kind(0.d0)),target :: work(:) - integer, intent(out) :: info - end subroutine psb_zmlprc_aply - end interface name='psb_zprc_aply' info = 0 @@ -180,7 +153,8 @@ end subroutine psb_zprc_aply !!$ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zprc_aply1 + implicit none type(psb_desc_type),intent(in) :: desc_data @@ -190,25 +164,12 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans) character(len=1), optional :: trans logical,parameter :: debug=.false., debugprt=.false. - interface - subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) - use psb_base_mod - use psb_prec_type - type(psb_desc_type),intent(in) :: desc_data - type(psb_zprec_type), intent(in) :: prec - complex(kind(0.d0)),intent(inout) :: x(:), y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(kind(0.d0)), optional, target :: work(:) - end subroutine psb_zprc_aply - end interface - ! Local variables character :: trans_ integer :: ictxt,np,me,i, isz, err_act, int_err(5) complex(kind(1.d0)), pointer :: WW(:), w1(:) character(len=20) :: name, ch_err - name='psb_zprec1' + name='psb_zprec_aply1' info = 0 call psb_erractionsave(err_act) diff --git a/psb_zprecbld.f90 b/psb_zprecbld.f90 index 7417838e..61201745 100644 --- a/psb_zprecbld.f90 +++ b/psb_zprecbld.f90 @@ -37,8 +37,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) use psb_base_mod - use psb_prec_type - use psb_prec_mod + use psb_prec_mod, mld_protect_name => psb_zprecbld Implicit None type(psb_zspmat_type), target :: a @@ -61,7 +60,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) info=0 err=0 call psb_erractionsave(err_act) - name = 'psb_precbld' + name = 'psb_zprecbld' if (debug) write(0,*) 'Entering precbld',P%prec,desc_a%matrix_data(:) info = 0 diff --git a/psb_zprecfree.f90 b/psb_zprecfree.f90 index 7dbe944a..84997aa9 100644 --- a/psb_zprecfree.f90 +++ b/psb_zprecfree.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_zprecfree(p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zprecfree + implicit none !....parameters... @@ -49,7 +50,7 @@ subroutine psb_zprecfree(p,info) if(psb_get_errstatus().ne.0) return info=0 - name = 'pszprecfree' + name = 'psb_zprecfree' call psb_erractionsave(err_act) me=-1 diff --git a/psb_zprecset.f90 b/psb_zprecset.f90 index 306d7352..a6d39b53 100644 --- a/psb_zprecset.f90 +++ b/psb_zprecset.f90 @@ -37,7 +37,8 @@ subroutine psb_zprecset(p,ptype,info,iv,rs,rv,ilev,nlev) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zprecset + implicit none type(psb_zprec_type), intent(inout) :: p diff --git a/psb_zslu_bld.f90 b/psb_zslu_bld.f90 index dd597402..140b5dfc 100644 --- a/psb_zslu_bld.f90 +++ b/psb_zslu_bld.f90 @@ -36,140 +36,46 @@ !!$ subroutine psb_zslu_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zslu_bld + implicit none - type(psb_zspmat_type), intent(inout) :: a + type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a type(psb_zbaseprc_type), intent(inout) :: p integer, intent(out) :: info - - type(psb_zspmat_type) :: blck, atmp - character(len=5) :: fmt - character :: upd='F' - integer :: i,j,nza,nzb,nzt,ictxt, me,np,err_act + integer :: i,j,nza,nzb,nzt,ictxt,me,np,err_act logical, parameter :: debug=.false. character(len=20) :: name, ch_err - interface psb_asmatbld - Subroutine psb_zasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_zspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_zasmatbld - end interface - if(psb_get_errstatus().ne.0) return info=0 - name='psb_slu_bld' + name='psb_zslu_bld' call psb_erractionsave(err_act) - ictxt = desc_A%matrix_data(psb_ctxt_) + ictxt = psb_cd_get_context(desc_a) call psb_info(ictxt, me, np) - fmt = 'COO' - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - atmp%fida='COO' - if (Debug) then - write(0,*) me, 'SPLUBLD: Calling csdp' - call psb_barrier(ictxt) - endif - - call psb_csdp(a,atmp,info) - if(info /= 0) then - info=4010 - ch_err='psb_csdp' - call psb_errpush(info,name,a_err=ch_err) + if (toupper(a%fida) /= 'CSR') then + write(0,*) 'Unimplemented input to SLU_BLD' goto 9999 - end if - nza = atmp%infoa(psb_nnz_) - if (Debug) then - write(0,*) me, 'SPLUBLD: Done csdp',info,nza,atmp%m,atmp%k - call psb_barrier(ictxt) endif - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=fmt) - if(info /= 0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - nzb = blck%infoa(psb_nnz_) - if (Debug) then - write(0,*) me, 'SPLUBLD: Done asmatbld',info,nzb,blck%fida - call psb_barrier(ictxt) - endif - if (nzb > 0 ) then - if (size(atmp%aspk) psb_zsp_renum + implicit none ! .. array Arguments .. @@ -57,7 +58,7 @@ subroutine psb_zsp_renum(a,desc_a,blck,p,atmp,info) if (psb_get_errstatus().ne.0) return info=0 - name='apply_renum' + name='psb_zsp_renum' call psb_erractionsave(err_act) ictxt=psb_cd_get_context(desc_a) diff --git a/psb_zumf_bld.f90 b/psb_zumf_bld.f90 index 9379dbb9..0a0be9d7 100644 --- a/psb_zumf_bld.f90 +++ b/psb_zumf_bld.f90 @@ -36,7 +36,8 @@ !!$ subroutine psb_zumf_bld(a,desc_a,p,info) use psb_base_mod - use psb_prec_type + use psb_prec_mod, mld_protect_name => psb_zumf_bld + implicit none type(psb_zspmat_type), intent(inout) :: a @@ -45,139 +46,43 @@ subroutine psb_zumf_bld(a,desc_a,p,info) integer, intent(out) :: info - type(psb_zspmat_type) :: blck, atmp - character(len=5) :: fmt - character :: upd='F' - integer :: i,j,nza,nzb,nzt,ictxt, me,np,err_act + integer :: i,j,nza,nzb,nzt,ictxt,me,np,err_act integer :: i_err(5) logical, parameter :: debug=.false. character(len=20) :: name, ch_err - interface psb_asmatbld - Subroutine psb_zasmatbld(ptype,novr,a,blk,desc_data,upd,desc_p,info,outfmt) - use psb_base_mod - use psb_prec_type - integer, intent(in) :: ptype,novr - Type(psb_zspmat_type), Intent(in) :: a - Type(psb_zspmat_type), Intent(inout) :: blk - Type(psb_desc_type), Intent(inout) :: desc_p - Type(psb_desc_type), Intent(in) :: desc_data - Character, Intent(in) :: upd - integer, intent(out) :: info - character(len=5), optional :: outfmt - end Subroutine psb_zasmatbld - end interface - info=0 - name='psb_umf_bld' + name='psb_zumf_bld' call psb_erractionsave(err_act) ictxt = desc_A%matrix_data(psb_ctxt_) call psb_info(ictxt, me, np) - fmt = 'COO' - call psb_nullify_sp(blck) - call psb_nullify_sp(atmp) - - atmp%fida='COO' - if (Debug) then - write(0,*) me, 'UMFBLD: Calling csdp' - call psb_barrier(ictxt) - endif - - call psb_zcsdp(a,atmp,info) - if(info /= 0) then - info=4010 - ch_err='psb_zcsdp' - call psb_errpush(info,name,a_err=ch_err) + if (toupper(a%fida) /= 'CSC') then + write(0,*) 'Unimplemented input to UMF_BLD' goto 9999 - end if - nza = psb_sp_get_nnzeros(atmp) - nzb = psb_sp_get_nnzeros(a) - if (Debug) then - write(0,*) me, 'UMFBLD: Done csdp',info,nza,atmp%m,atmp%k,nzb - call psb_barrier(ictxt) endif - call psb_asmatbld(p%iprcparm(p_type_),p%iprcparm(n_ovr_),a,& - & blck,desc_a,upd,p%desc_data,info,outfmt=fmt) - if(info /= 0) then - info=4010 - ch_err='psb_asmatbld' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - nzb = psb_sp_get_nnzeros(blck) - if (Debug) then - write(0,*) me, 'UMFBLD: Done asmatbld',info,nzb,blck%fida - call psb_barrier(ictxt) - endif - if (nzb > 0 ) then - if (size(atmp%aspk)