fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-20 12:52:14 +00:00
Code Issues Packages Projects Releases Wiki Activity

Completed reorganization of preconditioner build routines. Now it

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