fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Included first development version of multilevel stuff.

Browse Source
psblas3-type-indexed
Salvatore Filippone 18 years ago
parent 3f49de80d9
commit 048453938b
32 changed files with 1653 additions and 1001 deletions
Show Stats Download Patch File Download Diff File

7
src/modules/psb_prec_mod.f90
Unescape Escape View File

@ -39,7 +39,7 @@ module psb_prec_mod
use psb_prec_type
implicit none
type(psb_dspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_dprec_type), intent(inout) :: prec
integer, intent(out) :: info
character, intent(in),optional :: upd
@ -49,7 +49,7 @@ module psb_prec_mod
use psb_prec_type
implicit none
type(psb_zspmat_type), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_zprec_type), intent(inout) :: prec
integer, intent(out) :: info
character, intent(in),optional :: upd
@ -57,7 +57,7 @@ module psb_prec_mod
end interface
interface psb_precset
subroutine psb_dprecset(prec,ptype,info,iv,rs,rv)
subroutine psb_dprecset(prec,ptype,info,iv,rs,rv,ilev,nlev)
use psb_serial_mod
use psb_descriptor_type
use psb_prec_type
@ -66,6 +66,7 @@ module psb_prec_mod
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: iv(:)
integer, optional, intent(in) :: nlev,ilev
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
end subroutine psb_dprecset

134
src/modules/psb_prec_type.f90
Unescape Escape View File

@ -86,16 +86,60 @@ module psb_prec_type
type(psb_dspmat_type), pointer :: av(:) => null() !
real(kind(1.d0)), pointer :: d(:) => null()
type(psb_desc_type), pointer :: desc_data => null() !
type(psb_desc_type), pointer :: desc_data => null(), desc_ac=>null()! !
integer, pointer :: iprcparm(:) => null() !
real(kind(1.d0)), pointer :: dprcparm(:) => null() !
integer, pointer :: perm(:) => null(), invperm(:) => null()
integer, pointer :: mlia(:) => null(), nlaggr(:) => null() !
type(psb_dspmat_type), pointer :: aorig => null() !
type(psb_dspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc => null() !
real(kind(1.d0)), pointer :: dorig(:) => null() !
end type psb_dbaseprc_type
!
! Multilevel preconditioning
!
! To each level I there corresponds a matrix A(I) and a preconditioner K(I)
!
! A notational difference: in the DD reference above the preconditioner for
! a given level K(I) is written out as a sum over the subdomains
!
! SUM_k(R_k^T A_k R_k)
!
! whereas in this code the sum is implicit in the parallelization,
! i.e. each process takes care of one subdomain, and for each level we have
! as many subdomains as there are processes (except for the coarsest level where
! we might have a replicated index space). Thus the sum apparently disappears
! from our code, but only apparently, because it is implicit in the call
! to psb_baseprc_aply.
!
! A bit of description of the baseprecv(:) data structure:
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe
!
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev)
! if no smoother, it is used instead of sm_pr_
! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs.
!
type psb_dprec_type
type(psb_dbaseprc_type), pointer :: baseprecv(:) => null()
! contain type of preconditioning to be performed
@ -106,12 +150,13 @@ module psb_prec_type
type(psb_zspmat_type), pointer :: av(:) => null() !
complex(kind(1.d0)), pointer :: d(:) => null()
type(psb_desc_type), pointer :: desc_data => null() !
type(psb_desc_type), pointer :: desc_data => null() , desc_ac=>null()!
integer, pointer :: iprcparm(:) => null() !
real(kind(1.d0)), pointer :: dprcparm(:) => null() !
integer, pointer :: perm(:) => null(), invperm(:) => null()
integer, pointer :: mlia(:) => null(), nlaggr(:) => null() !
type(psb_zspmat_type), pointer :: aorig => null() !
type(psb_zspmat_type), pointer :: base_a => null() !
type(psb_desc_type), pointer :: base_desc => null() !
complex(kind(1.d0)), pointer :: dorig(:) => null() !
end type psb_zbaseprc_type
@ -181,6 +226,7 @@ contains
subroutine psb_file_prec_descr(iout,p)
integer, intent(in) :: iout
type(psb_dprec_type), intent(in) :: p
integer :: ilev
write(iout,*) 'Preconditioner description'
if (associated(p%baseprecv)) then
@ -206,40 +252,44 @@ contains
end select
end if
if (size(p%baseprecv)>=2) then
if (.not.associated(p%baseprecv(2)%iprcparm)) then
do ilev = 2, size(p%baseprecv)
if (.not.associated(p%baseprecv(ilev)%iprcparm)) then
write(iout,*) 'Inconsistent MLPREC part!'
return
endif
write(iout,*) 'Multilevel: ',ml_names(p%baseprecv(2)%iprcparm(ml_type_))
if (p%baseprecv(2)%iprcparm(ml_type_)>no_ml_) then
write(iout,*) 'Multilevel: Level No', ilev
write(iout,*) 'Multilevel type: ',&
& ml_names(p%baseprecv(ilev)%iprcparm(ml_type_))
if (p%baseprecv(ilev)%iprcparm(ml_type_)>no_ml_) then
write(iout,*) 'Multilevel aggregation: ', &
& aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
& aggr_names(p%baseprecv(ilev)%iprcparm(aggr_alg_))
write(iout,*) 'Smoother: ', &
& smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_))
if (p%baseprecv(2)%iprcparm(smth_kind_) /= no_smth_) then
write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
& smooth_kinds(p%baseprecv(ilev)%iprcparm(smth_kind_))
if (p%baseprecv(ilev)%iprcparm(smth_kind_) /= no_smth_) then
write(iout,*) 'Smoothing omega: ', p%baseprecv(ilev)%dprcparm(smooth_omega_)
write(iout,*) 'Smoothing position: ',&
& smooth_names(p%baseprecv(2)%iprcparm(smth_pos_))
& smooth_names(p%baseprecv(ilev)%iprcparm(smth_pos_))
end if
write(iout,*) 'Coarse matrix: ',&
& matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
& matrix_names(p%baseprecv(ilev)%iprcparm(coarse_mat_))
write(iout,*) 'Aggregation sizes: ', &
& sum( p%baseprecv(2)%nlaggr(:)),' : ',p%baseprecv(2)%nlaggr(:)
& sum( p%baseprecv(ilev)%nlaggr(:)),' : ',p%baseprecv(ilev)%nlaggr(:)
write(iout,*) 'Factorization type: ',&
& fact_names(p%baseprecv(2)%iprcparm(f_type_))
select case(p%baseprecv(2)%iprcparm(f_type_))
& fact_names(p%baseprecv(ilev)%iprcparm(f_type_))
select case(p%baseprecv(ilev)%iprcparm(f_type_))
case(f_ilu_n_)
write(iout,*) 'Fill level :',p%baseprecv(2)%iprcparm(ilu_fill_in_)
write(iout,*) 'Fill level :',p%baseprecv(ilev)%iprcparm(ilu_fill_in_)
case(f_ilu_e_)
write(iout,*) 'Fill threshold :',p%baseprecv(2)%dprcparm(fact_eps_)
write(iout,*) 'Fill threshold :',p%baseprecv(ilev)%dprcparm(fact_eps_)
case(f_slu_,f_umf_)
case default
write(iout,*) 'Should never get here!'
end select
write(iout,*) 'Number of Jacobi sweeps: ', &
& (p%baseprecv(2)%iprcparm(jac_sweeps_))
& (p%baseprecv(ilev)%iprcparm(jac_sweeps_))
end if
end do
end if
else
@ -355,13 +405,14 @@ contains
& aggr_names(p%baseprecv(2)%iprcparm(aggr_alg_))
write(iout,*) 'Smoother: ', &
& smooth_kinds(p%baseprecv(2)%iprcparm(smth_kind_))
write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
if (p%baseprecv(2)%iprcparm(smth_kind_) /= no_smth_) then
write(iout,*) 'Smoothing omega: ', p%baseprecv(2)%dprcparm(smooth_omega_)
write(iout,*) 'Smoothing position: ',&
& smooth_names(p%baseprecv(2)%iprcparm(smth_pos_))
end if
write(iout,*) 'Coarse matrix: ',&
& matrix_names(p%baseprecv(2)%iprcparm(coarse_mat_))
end if
write(iout,*) 'Aggregation sizes: ', &
& sum( p%baseprecv(2)%nlaggr(:)),' : ',p%baseprecv(2)%nlaggr(:)
write(iout,*) 'Factorization type: ',&
@ -631,12 +682,23 @@ contains
end if
deallocate(p%desc_data)
endif
if (associated(p%desc_ac)) then
if (associated(p%desc_ac%matrix_data)) then
call psb_cdfree(p%desc_ac,info)
end if
deallocate(p%desc_ac)
endif
if (associated(p%dprcparm)) then
deallocate(p%dprcparm,stat=info)
end if
if (associated(p%aorig)) then
if (associated(p%base_a)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
endif
if (associated(p%base_desc)) then
! This is a pointer to something else, must not free it here.
nullify(p%aorig)
nullify(p%base_desc)
endif
if (associated(p%dorig)) then
deallocate(p%dorig,stat=info)
@ -676,7 +738,7 @@ contains
type(psb_dbaseprc_type), intent(inout) :: p
nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
& p%nlaggr,p%aorig,p%dorig,p%desc_data)
& p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data, p%desc_ac)
end subroutine psb_nullify_dbaseprec
@ -712,12 +774,22 @@ contains
end if
deallocate(p%desc_data)
endif
if (associated(p%desc_ac)) then
if (associated(p%desc_ac%matrix_data)) then
call psb_cdfree(p%desc_ac,info)
end if
deallocate(p%desc_ac)
endif
if (associated(p%dprcparm)) then
deallocate(p%dprcparm,stat=info)
end if
if (associated(p%aorig)) then
if (associated(p%base_a)) then
! This is a pointer to something else, must not free it here.
nullify(p%base_a)
endif
if (associated(p%base_desc)) then
! This is a pointer to something else, must not free it here.
nullify(p%aorig)
nullify(p%base_desc)
endif
if (associated(p%dorig)) then
deallocate(p%dorig,stat=info)
@ -741,11 +813,11 @@ contains
if (associated(p%iprcparm)) then
if (p%iprcparm(f_type_)==f_slu_) then
!!$ call psb_zslu_free(p%iprcparm(slu_ptr_),info)
call psb_zslu_free(p%iprcparm(slu_ptr_),info)
end if
if (p%iprcparm(f_type_)==f_umf_) then
!!$ call psb_zumf_free(p%iprcparm(umf_symptr_),&
!!$ & p%iprcparm(umf_numptr_),info)
call psb_zumf_free(p%iprcparm(umf_symptr_),&
& p%iprcparm(umf_numptr_),info)
end if
deallocate(p%iprcparm,stat=info)
end if
@ -757,7 +829,7 @@ contains
type(psb_zbaseprc_type), intent(inout) :: p
nullify(p%av,p%d,p%iprcparm,p%dprcparm,p%perm,p%invperm,p%mlia,&
& p%nlaggr,p%aorig,p%dorig,p%desc_data)
& p%nlaggr,p%base_a,p%base_desc,p%dorig,p%desc_data,p%desc_ac)
end subroutine psb_nullify_zbaseprec

6
src/modules/psb_serial_mod.f90
Unescape Escape View File

@ -359,13 +359,15 @@ module psb_serial_mod
interface psb_symbmm
subroutine psb_dsymbmm(a,b,c)
subroutine psb_dsymbmm(a,b,c,info)
use psb_spmat_type
type(psb_dspmat_type) :: a,b,c
integer :: info
end subroutine psb_dsymbmm
subroutine psb_zsymbmm(a,b,c)
subroutine psb_zsymbmm(a,b,c,info)
use psb_spmat_type
type(psb_zspmat_type) :: a,b,c
integer :: info
end subroutine psb_zsymbmm
end interface

67
src/prec/psb_dbaseprc_aply.f90
Unescape Escape View File

@ -34,9 +34,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
@ -52,7 +52,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type),intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: prec
real(kind(0.d0)),intent(inout) :: x(:), y(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -68,13 +68,13 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name, ch_err
interface psb_bjac_aply
subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type), intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: prec
real(kind(0.d0)),intent(inout) :: x(:), y(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -105,33 +105,35 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
case(noprec_)
n_row=desc_data%matrix_data(psb_n_row_)
if (beta==dzero) then
y(1:n_row) = x(1:n_row)
else if (beta==done) then
y(1:n_row) = x(1:n_row) + y(1:n_row)
else if (beta==-done) then
y(1:n_row) = x(1:n_row) - y(1:n_row)
else
y(1:n_row) = x(1:n_row) + beta*y(1:n_row)
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(diagsc_)
n_row=desc_data%matrix_data(psb_n_row_)
if (beta==dzero) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row)
else if (beta==done) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row)
else if (beta==-done) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row)
if (size(work) >= size(x)) then
ww => work
else
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row)
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
end if
n_row=desc_data%matrix_data(psb_n_row_)
ww(1:n_row) = x(1:n_row)*prec%d(1:n_row)
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (size(work) < size(x)) then
deallocate(ww,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Deallocate')
goto 9999
end if
end if
case(bja_)
call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjac_aply'
@ -142,7 +144,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
if (prec%iprcparm(n_ovr_)==0) then
! shortcut: this fixes performance for RAS(0) == BJA
call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjacaply'
@ -214,7 +216,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
end if
endif
call psb_bjac_aply(prec,tx,dzero,ty,prec%desc_data,trans,aux,info)
call psb_bjac_aply(done,prec,tx,dzero,ty,prec%desc_data,trans,aux,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjac_aply'
@ -250,18 +252,7 @@ subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
& prec%iprcparm(prol_)
end select
if (beta == dzero) then
y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_))
else if (beta == done) then
y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
else if (beta == -done) then
y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
else
y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
if ((6*isz) <= size(work)) then

13
src/prec/psb_dbaseprc_bld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd)
Implicit None
type(psb_dspmat_type), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_dbaseprc_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -169,7 +169,13 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd)
select case(p%iprcparm(p_type_))
case (noprec_)
! Do nothing.
call psb_cdcpy(desc_a,p%desc_data,info)
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (diagsc_)
@ -256,7 +262,8 @@ subroutine psb_dbaseprc_bld(a,desc_a,p,info,upd)
end select
p%base_a => a
p%base_desc => desc_a
call psb_erractionrestore(err_act)
return

41
src/prec/psb_dbjac_aply.f90
Unescape Escape View File

@ -34,9 +34,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a Block Jacobi preconditioner stored in prec
! Note that desc_data may or may not be the same as prec%desc_data,
! but since both are INTENT(IN) this should be legal.
@ -54,7 +54,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type), intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: prec
real(kind(0.d0)),intent(inout) :: x(:), y(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -125,7 +125,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
& trans='N',unit=diagl,choice=psb_none_,work=aux)
if(info /=0) goto 9999
ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row)
call psb_spsm(done,prec%av(u_pr_),ww,beta,y,desc_data,info,&
call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit=diagu,choice=psb_none_, work=aux)
if(info /=0) goto 9999
@ -134,7 +134,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
& trans=trans,unit=diagu,choice=psb_none_, work=aux)
if(info /=0) goto 9999
ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row)
call psb_spsm(done,prec%av(l_pr_),ww,beta,y,desc_data,info,&
call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit=diagl,choice=psb_none_,work=aux)
if(info /=0) goto 9999
@ -152,16 +152,8 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (beta == dzero) then
y(1:n_row) = ww(1:n_row)
else if (beta==done) then
y(1:n_row) = ww(1:n_row) + y(1:n_row)
else if (beta==-done) then
y(1:n_row) = ww(1:n_row) - y(1:n_row)
else
y(1:n_row) = ww(1:n_row) + beta*y(1:n_row)
endif
case (f_umf_)
@ -174,15 +166,7 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
if (beta == dzero) then
y(1:n_row) = ww(1:n_row)
else if (beta==dzero) then
y(1:n_row) = ww(1:n_row) + y(1:n_row)
else if (beta==-dzero) then
y(1:n_row) = ww(1:n_row) - y(1:n_row)
else
y(1:n_row) = ww(1:n_row) + beta*y(1:n_row)
endif
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_)
@ -253,15 +237,8 @@ subroutine psb_dbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
end select
if (beta == dzero) then
y(1:n_row) = tx(1:n_row)
else if (beta==done) then
y(1:n_row) = tx(1:n_row) + y(1:n_row)
else if (beta==-done) then
y(1:n_row) = tx(1:n_row) - y(1:n_row)
else
y(1:n_row) = tx(1:n_row) + beta*y(1:n_row)
endif
call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
deallocate(tx,ty)

41
src/prec/psb_dbldaggrmat.f90
Unescape Escape View File

@ -75,13 +75,14 @@ subroutine psb_dbldaggrmat(a,desc_a,ac,p,desc_p,info)
if (aggr_dump) call psb_csprt(90+me,ac,head='% Raw aggregate.')
case(smth_omg_,smth_biz_)
if (aggr_dump) call psb_csprt(70+me,a,head='% Input matrix')
call smooth_aggregate(info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='smooth_aggregate')
goto 9999
end if
if (aggr_dump) call psb_csprt(90+me,ac,head='% Smooth aggregate.')
case default
call psb_errpush(4010,name,a_err=name)
goto 9999
@ -117,6 +118,7 @@ contains
integer :: ictxt, nrow, nglob, ncol, ntaggr, nzbg, ip, ndx,&
& naggr, np, me, nzt,irs,jl,nzl,nlr,&
& icomm,naggrm1, i, j, k, err_act
name='raw_aggregate'
if(psb_get_errstatus().ne.0) return
info=0
@ -228,11 +230,6 @@ contains
enddo
end if
call psb_fixcoo(b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='fixcoo')
goto 9999
end if
irs = b%infoa(psb_nnz_)
call psb_sp_reall(b,irs,info)
@ -523,8 +520,6 @@ contains
endif
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
call psb_ipcoo2csr(am4,info)
@ -542,7 +537,7 @@ contains
!
! WARNING: the cycles below assume that AM3 does have
! its diagonal elements stored explicitly!!!
! Should we swicth to something safer?
! Should we switch to something safer?
!
call psb_spscal(am3,p%dorig,info)
if(info /= 0) goto 9999
@ -604,12 +599,15 @@ contains
am3%aspk(j) = done - omega*am3%aspk(j)
endif
end do
call psb_ipcoo2csr(am3,info)
else
write(0,*) 'Missing implementation of I sum'
call psb_errpush(4010,name)
goto 9999
end if
if (test_dump) call &
& psb_csprt(20+me,am4,head='% Operator Ptilde.',ivr=desc_a%loc_to_glob)
if (test_dump) call psb_csprt(40+me,am3,head='% (I-wDA)',ivr=desc_a%loc_to_glob,&
& ivc=desc_a%loc_to_glob)
if (debug) write(0,*) me,'Done gather, going for SYMBMM 1'
@ -620,7 +618,15 @@ contains
! Doing it this way means to consider diag(Ai)
!
!
call psb_symbmm(am3,am4,am1)
call psb_symbmm(am3,am4,am1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 1')
goto 9999
end if
am1%aspk(:) = 0.d0
if (test_dump) &
& call psb_csprt(50+me,am1,head='% (I-wDA)Pt symbmm ')
call psb_numbmm(am3,am4,am1)
if (debug) write(0,*) me,'Done NUMBMM 1'
@ -667,7 +673,12 @@ contains
if (test_dump) &
& call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob)
call psb_symbmm(a,am1,am3)
call psb_symbmm(a,am1,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 2')
goto 9999
end if
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
@ -724,7 +735,12 @@ contains
endif
if (debug) write(0,*) me,'starting symbmm 3'
call psb_symbmm(am2,am3,b)
call psb_symbmm(am2,am3,b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm 3')
goto 9999
end if
if (debug) write(0,*) me,'starting numbmm 3'
call psb_numbmm(am2,am3,b)
if (debug) write(0,*) me,'Done NUMBMM 3'
@ -1018,6 +1034,7 @@ contains
deallocate(nzbr,idisp)
end select
call psb_ipcoo2csr(bg,info)
call psb_erractionrestore(err_act)
return

5
src/prec/psb_ddiagsc_bld.f90
Unescape Escape View File

@ -98,6 +98,11 @@ subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psb_cdcpy(desc_a,p%desc_Data,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdcpy')
goto 9999
end if
if (debug) write(ilout+me,*) 'VDIAG ',n_row
do i=1,n_row

692
src/prec/psb_dmlprc_aply.f90
Unescape Escape View File

@ -34,10 +34,53 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! where K is a multilevel (actually 2-level) preconditioner stored in prec
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a multilevel preconditioner stored in baseprecv
!
! cfr.: Smith, Biorstad & Gropp
! Domain Decomposition
! Cambridge Univ. Press
!
! To each level I there corresponds a matrix A(I) and a preconditioner K(I)
!
! A notational difference: in the DD reference above the preconditioner for
! a given level K(I) is written out as a sum over the subdomains
!
! SUM_k(R_k^T A_k R_k)
!
! whereas in this code the sum is implicit in the parallelization,
! i.e. each process takes care of one subdomain, and for each level we have
! as many subdomains as there are processes (except for the coarsest level where
! we might have a replicated index space). Thus the sum apparently disappears
! from our code, but only apparently, because it is implicit in the call
! to psb_baseprc_aply.
!
! A bit of description of the baseprecv(:) data structure:
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe
!
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev)
! if no smoother, it is used instead of sm_pr_
! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs.
!
use psb_serial_mod
@ -52,41 +95,49 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type),intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: baseprecv(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
real(kind(0.d0)),intent(inout) :: x(:), y(:)
character :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,n_col
real(kind(1.d0)), allocatable :: tx(:),ty(:),t2l(:),w2l(:),&
& x2l(:),b2l(:),tz(:),tty(:)
real(kind(1.d0)), allocatable :: tx(:),ty(:),x2l(:),y2l(:),&
& b2l(:),tty(:)
character ::diagl, diagu
integer :: ictxt,np,me,i, isz, nrg,nr2l,err_act, iptype, int_err(5)
real(kind(1.d0)) :: omega
real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime
logical, parameter :: debug=.false., debugprt=.false.
integer :: ismth
integer :: ismth, nlev, ilev
external mpi_wtime
character(len=20) :: name, ch_err
type psb_mlprec_wrk_type
real(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),&
& x2l(:)=>null(),y2l(:)=>null(),&
& b2l(:)=>null(),tty(:)=>null()
end type psb_mlprec_wrk_type
type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:)
interface psb_baseprc_aply
subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: prec
real(kind(0.d0)),intent(inout) :: x(:), y(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
end subroutine psb_dbaseprc_aply
end interface
name='psb_dmlprc_aply'
name='psb_mlprc_aply'
info = 0
call psb_erractionsave(err_act)
@ -94,426 +145,625 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
omega=baseprecv(2)%dprcparm(smooth_omega_)
ismth=baseprecv(2)%iprcparm(smth_kind_)
nlev = size(baseprecv)
allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
select case(baseprecv(2)%iprcparm(ml_type_))
case(no_ml_)
! Should not really get here.
write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????'
call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?')
goto 9999
case(add_ml_prec_)
!
! Additive multilevel
! Additive is very simple.
! 1. X(1) = Xext
! 2. DO ILEV=2,NLEV
! X(ILEV) = AV(PR_SM_T_)*X(ILEV-1)
! 3. Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! 4. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = AV(PR_SM_)*Y(ILEV+1)
! 5. Yext = beta*Yext + Y(1)
!
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),stat=info)
call psb_baseprc_aply(alpha,baseprecv(1),x,beta,y,&
& baseprecv(1)%base_desc,trans,work,info)
if(info /=0) goto 9999
allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y)))
mlprec_wrk(1)%x2l(:) = x(:)
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
allocate(mlprec_wrk(ilev)%x2l(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%tx(max(n_row,n_col)),&
& mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = dzero
w2l(:) = dzero
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(1:n_row) = mlprec_wrk(ilev-1)%x2l(1:n_row)
mlprec_wrk(ilev)%tx(n_row+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev)%ty(:) = dzero
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
!
! Smoothed aggregation
!
allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),&
& tz(max(n_row,n_col)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_))
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = dzero
end if
call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info)
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcut
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i)
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%x2l(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_baseprc_aply(done,baseprecv(ilev),&
& mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev)%desc_data, 'N',work,info)
enddo
do ilev =nlev,2,-1
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
if (ismth /= no_smth_) then
call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info)
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,info)
if(info /=0) goto 9999
!
! Finally add back into Y.
!
call psb_geaxpby(done,ty,done,y,desc_data,info)
if(info /=0) goto 9999
deallocate(tx,ty,tz)
else
do i=1, desc_data%matrix_data(psb_n_row_)
y(i) = y(i) + t2l(baseprecv(2)%mlia(i))
do i=1, n_row
mlprec_wrk(ilev-1)%y2l(i) = mlprec_wrk(ilev-1)%y2l(i) + &
& mlprec_wrk(ilev)%y2l(baseprecv(ilev)%mlia(i))
enddo
end if
end do
if (debugprt) write(0,*)' Y2: ',Y(:)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,done,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
deallocate(t2l,w2l)
case(mult_ml_prec_)
!
! Multiplicative multilevel
! Pre/post smoothing versions.
!
select case(baseprecv(2)%iprcparm(smth_pos_))
case(post_smooth_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = dzero
w2l(:) = dzero
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==dzero
! Post smoothing.
! 1. X(1) = Xext
! 2. DO ILEV=2, NLEV :: X(ILEV) = AV(PR_SM_T_,ILEV)*X(ILEV-1)
! 3. Y(NLEV) = (K(NLEV)**(-1))*X(NLEV)
! 4. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV))
!
! 5. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
! Also: post smoothing is not spelled out in detail in DD.
!
!
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%tx(nr2l), stat=info)
mlprec_wrk(1)%x2l(:) = dzero
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%tx(:) = dzero
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
& baseprecv(1)%base_desc,info)
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
do ilev=2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth = baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
if (debug) write(0,*)' mult_ml_apply omega ',omega
if (debugprt) write(0,*)' mult_ml_apply X: ',X(:)
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
if (info /= 0) then
if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg
call psb_errpush(4010,name,a_err='axpby post_smooth 1')
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= no_smth_) then
!
! Smoothed aggregation
!
allocate(tz(max(n_row,n_col)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = dzero
end if
call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info)
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& dzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcut
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%x2l(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
enddo
call psb_baseprc_aply(done,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
& dzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,&
&'N',work,info)
if(info /=0) goto 9999
do ilev=nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info)
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& dzero,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
!
! Finally add back into Y.
!
deallocate(tz)
else
ty(:) = dzero
do i=1, desc_data%matrix_data(psb_n_row_)
ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i))
n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_)
mlprec_wrk(ilev)%y2l(:) = dzero
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
deallocate(t2l,w2l)
call psb_spmm(-done,baseprecv(2)%aorig,ty,done,tx,desc_data,info,work=work)
if(info /=0) goto 9999
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
call psb_baseprc_aply(baseprecv(1),tx,done,ty,desc_data,trans,&
& work,info)
if(info /=0) goto 9999
call psb_geaxpby(done,ty,beta,y,desc_data,info)
call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
deallocate(tx,ty)
enddo
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
case(pre_smooth_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info)
!
! Pre smoothing.
! 1. X(1) = Xext
! 2. Y(1) = (K(1)**(-1))*X(1)
! 3. TX(1) = X(1) - A(1)*Y(1)
! 4. DO ILEV=2, NLEV
! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1)
! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV))
! 5. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! 6. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
!
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%tx(nr2l), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = dzero
w2l(:) = dzero
mlprec_wrk(1)%y2l(:) = dzero
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==zero
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if(info /=0) goto 9999
call psb_baseprc_aply(baseprecv(1),x,dzero,tty,desc_data,&
mlprec_wrk(1)%x2l(:) = x
call psb_baseprc_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,&
& dzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work)
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (ismth /= no_smth_) then
allocate(tz(max(n_row,n_col)),stat=info)
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth = baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (ismth /= no_smth_) then
!
!Smoothed Aggregation
!
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev-1)%tx(n_row+1:max(n_row,n_col)) = dzero
end if
call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info)
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcuts
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
mlprec_wrk(ilev)%x2l = dzero
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%tx(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /=0) goto 9999
if(ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
endif
enddo
do ilev = nlev-1, 1, -1
ismth=baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info)
if(info /=0) goto 9999
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,info)
call psb_geaxpby(done,ty,done,tty,desc_data,info)
if(info /=0) goto 9999
deallocate(tz)
else
do i=1, desc_data%matrix_data(psb_n_row_)
tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i))
n_row = baseprecv(ilev+1)%base_desc%matrix_data(psb_n_row_)
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
call psb_geaxpby(done,tty,beta,y,desc_data,info)
enddo
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
deallocate(t2l,w2l,tx,ty,tty)
case(smooth_both_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info)
!
! Symmetrized smoothing.
! 1. X(1) = Xext
! 2. Y(1) = (K(1)**(-1))*X(1)
! 3. TX(1) = X(1) - A(1)*Y(1)
! 4. DO ILEV=2, NLEV
! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1)
! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV))
! 5. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV))
! 6. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
!
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%ty(nr2l), mlprec_wrk(1)%tx(nr2l), stat=info)
mlprec_wrk(1)%x2l(:) = dzero
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%tx(:) = dzero
mlprec_wrk(1)%ty(:) = dzero
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = dzero
w2l(:) = dzero
tx(:) = dzero
ty(:) = dzero
tty(:) = dzero
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%tx,&
& baseprecv(1)%base_desc,info)
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==zero
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if(info /=0) goto 9999
call psb_baseprc_aply(done,baseprecv(1),mlprec_wrk(1)%x2l,&
& dzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
call psb_baseprc_aply(baseprecv(1),tx,dzero,tty,desc_data,trans,work,info)
if(info /=0) goto 9999
call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work)
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
call psb_spmm(-done,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& done,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%ty(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
mlprec_wrk(ilev)%ty(:) = dzero
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
if (ismth /= no_smth_) then
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,baseprecv(1)%desc_data,info,work=work)
!
!Smoothed Aggregation
!
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = dzero
mlprec_wrk(ilev-1)%ty(n_row+1:max(n_row,n_col)) = dzero
end if
call psb_csmm(done,baseprecv(2)%av(sm_pr_t_),tx,dzero,t2l,info)
call psb_csmm(done,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,dzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcuts
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
mlprec_wrk(ilev)%x2l = dzero
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%ty(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,dzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_geaxpby(done,mlprec_wrk(ilev)%x2l,dzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
if (ismth /= no_smth_) then
call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& dzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(done,baseprecv(2)%av(sm_pr_),t2l,dzero,ty,info)
if(info /=0) goto 9999
call psb_geaxpby(done,ty,done,tty,desc_data,info)
if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
endif
else
enddo
do ilev=nlev-1, 1, -1
ismth=baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(done,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
do i=1, desc_data%matrix_data(psb_n_row_)
tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i))
else
n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_)
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_spmm(-done,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
call psb_spmm(-done,baseprecv(2)%aorig,tty,done,tx,desc_data,info,work=work)
call psb_baseprc_aply(done,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& done,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
call psb_baseprc_aply(baseprecv(1),tx,done,tty,desc_data,'N',work,info)
enddo
call psb_geaxpby(done,tty,beta,y,desc_data,info)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
deallocate(t2l,w2l,tx,ty,tty)
if(info /=0) goto 9999
case default
write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_)
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(smth_pos_),0,0,0,0/))
goto 9999
end select
case default
write(0,*) me, 'Wrong mltype into PRC_APLY ',&
& baseprecv(2)%iprcparm(ml_type_)
call psb_errpush(4013,name,a_err='wrong mltype',&
& i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/))
goto 9999
end select
call mlprec_wrk_free(mlprec_wrk)
deallocate(mlprec_wrk)
call psb_erractionrestore(err_act)
return
@ -526,4 +776,24 @@ subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
end if
return
contains
subroutine mlprec_wrk_free(wrk)
type(psb_mlprec_wrk_type) :: wrk(:)
! This will not be needed when we have allocatables, as
! it is sufficient to deallocate the container, and
! the compiler is supposed to recursively deallocate the
! various components.
integer i
do i=1, size(wrk)
if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
end do
end subroutine mlprec_wrk_free
end subroutine psb_dmlprc_aply

14
src/prec/psb_dmlprc_bld.f90
Unescape Escape View File

@ -135,8 +135,6 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
p%aorig => a
nullify(p%d)
@ -167,7 +165,7 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
call psb_baseprc_bld(ac,desc_p,p,info)
if (debug) write(0,*) 'Out from basaeprcbld',info
if (debug) write(0,*) 'Out from baseprcbld',info
if(info /= 0) then
info=4010
ch_err='psb_baseprc_bld'
@ -180,13 +178,13 @@ subroutine psb_dmlprc_bld(a,desc_a,p,info)
! We have used a separate ac because:
! 1. We want to reuse the same routines psb_ilu_bld etc.
! 2. We do NOT want to pass an argument twice to them
! p%av(ac_) and p
! Hence a separate AC and a TRANSFER function.
! p%av(ac_) and p, as this would violate the Fortran standard
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)
call psb_cdfree(desc_p,info)
deallocate(desc_p)
p%base_a => p%av(ac_)
p%desc_ac => desc_p
nullify(desc_p)
call psb_erractionrestore(err_act)
return

12
src/prec/psb_dprc_aply.f90
Unescape Escape View File

@ -61,13 +61,13 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
character(len=20) :: name
interface psb_baseprc_aply
subroutine psb_dbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_dbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: prec
real(kind(0.d0)),intent(inout) :: x(:), y(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
real(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -75,12 +75,12 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end interface
interface psb_mlprc_aply
subroutine psb_dmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
subroutine psb_dmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_dbaseprc_type), intent(in) :: baseprecv(:)
real(kind(0.d0)),intent(in) :: beta
real(kind(0.d0)),intent(in) :: alpha,beta
real(kind(0.d0)),intent(inout) :: x(:), y(:)
character :: trans
real(kind(0.d0)),target :: work(:)
@ -117,14 +117,14 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end if
if (size(prec%baseprecv) >1) then
if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y)
call psb_mlprc_aply(prec%baseprecv,x,dzero,y,desc_data,trans_,work_,info)
call psb_mlprc_aply(done,prec%baseprecv,x,dzero,y,desc_data,trans_,work_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_dmlprc_aply')
goto 9999
end if
else if (size(prec%baseprecv) == 1) then
call psb_baseprc_aply(prec%baseprecv(1),x,dzero,y,desc_data,trans_, work_,info)
call psb_baseprc_aply(done,prec%baseprecv(1),x,dzero,y,desc_data,trans_, work_,info)
else
write(0,*) 'Inconsistent preconditioner: size of baseprecv???'
endif

36
src/prec/psb_dprecbld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
Implicit None
type(psb_dspmat_type), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_dprec_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -59,8 +59,8 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
Use psb_spmat_type
use psb_descriptor_type
use psb_prec_type
type(psb_dspmat_type) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), target :: a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_dbaseprc_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -149,7 +149,9 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
endif
if (size(p%baseprecv) > 1) then
call init_baseprc_av(p%baseprecv(2),info)
do i=2, size(p%baseprecv)
call init_baseprc_av(p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='allocate'
@ -157,31 +159,17 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
goto 9999
endif
call psb_mlprc_bld(a,desc_a,p%baseprecv(2),info)
call psb_mlprc_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,&
& p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='psb_mlprc_bld'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name)
goto 9999
endif
!
! 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
if (debug) then
write(0,*) 'Return from ',i-1,' call to mlprcbld ',info
endif
call psb_mlprc_bld(p%baseprecv(i-1)%av(ac_),p%baseprecv(i-1)%desc_data,&
& p%baseprecv(i),info)
end do
endif
@ -206,7 +194,7 @@ contains
! Have not decided what to do yet
end if
allocate(p%av(max_avsz),stat=info)
if (info /= 0) return
!!$ if (info /= 0) return
do k=1,size(p%av)
call psb_nullify_sp(p%av(k))
end do

209
src/prec/psb_dprecset.f90
Unescape Escape View File

@ -34,33 +34,57 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dprecset(p,ptype,info,iv,rs,rv)
subroutine psb_dprecset(p,ptype,info,iv,rs,rv,ilev,nlev)
use psb_serial_mod
use psb_descriptor_type
use psb_prec_type
implicit none
type(psb_dprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype
integer, intent(out) :: info
integer, optional, intent(in) :: iv(:)
integer, optional, intent(in) :: nlev,ilev
real(kind(1.d0)), optional, intent(in) :: rs
real(kind(1.d0)), optional, intent(in) :: rv(:)
type(psb_dbaseprc_type), pointer :: bpv(:)=>null()
character(len=len(ptype)) :: typeup
integer :: isz, err
integer :: isz, err, nlev_, ilev_, i
info = 0
if (present(ilev)) then
ilev_ = max(1, ilev)
else
ilev_ = 1
end if
if (present(nlev)) then
if (associated(p%baseprecv)) then
write(0,*) 'Warning: NLEV is ignored when P is already associated'
end if
nlev_ = max(1, nlev)
else
nlev_ = 1
end if
if (.not.associated(p%baseprecv)) then
allocate(p%baseprecv(1),stat=err)
call psb_nullify_baseprec(p%baseprecv(1))
allocate(p%baseprecv(nlev_),stat=err)
do i=1, nlev_
call psb_nullify_baseprec(p%baseprecv(i))
end do
else
nlev_ = size(p%baseprecv)
endif
if ((ilev_<1).or.(ilev_ > nlev_)) then
write(0,*) 'PRECSET ERRROR: ilev out of bounds'
info = -1
return
endif
if (.not.associated(p%baseprecv(1)%iprcparm)) then
allocate(p%baseprecv(1)%iprcparm(ifpsz),stat=err)
if (.not.associated(p%baseprecv(ilev_)%iprcparm)) then
allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),&
& p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure',err
endif
@ -68,124 +92,103 @@ subroutine psb_dprecset(p,ptype,info,iv,rs,rv)
select case(toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = noprec_
p%baseprecv(1)%iprcparm(f_type_) = f_none_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = noprec_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('DIAG','DIAGSC')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = diagsc_
p%baseprecv(1)%iprcparm(f_type_) = f_none_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = diagsc_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_none_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('BJA','ILU')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(1)%iprcparm(p_type_) = bja_
p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(1)%iprcparm(restr_) = psb_none_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 0
p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bja_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
case ('ASM','AS')
p%baseprecv(1)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(:) = 0
! Defaults first
p%baseprecv(1)%iprcparm(p_type_) = asm_
p%baseprecv(1)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(1)%iprcparm(restr_) = psb_halo_
p%baseprecv(1)%iprcparm(prol_) = psb_none_
p%baseprecv(1)%iprcparm(iren_) = 0
p%baseprecv(1)%iprcparm(n_ovr_) = 1
p%baseprecv(1)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(1)%iprcparm(jac_sweeps_) = 1
p%baseprecv(ilev_)%iprcparm(p_type_) = asm_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_halo_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 1
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(1)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(1)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(1)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(1)%iprcparm(f_type_) = iv(4)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(n_ovr_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(restr_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(prol_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(4)
! Do not consider renum for the time being.
!!$ if (isz >= 5) p%baseprecv(1)%iprcparm(iren_) = iv(5)
!!$ if (isz >= 5) p%baseprecv(ilev_)%iprcparm(iren_) = iv(5)
end if
case ('ML', '2L', '2LEV')
select case (size(p%baseprecv))
case(1)
! Reallocate
allocate(bpv(2),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:1',err
endif
bpv(1) = p%baseprecv(1)
call psb_nullify_baseprec(bpv(2))
deallocate(p%baseprecv)
p%baseprecv => bpv
nullify(bpv)
case(2)
! Do nothing
case default
! Error
end select
allocate(p%baseprecv(2)%iprcparm(ifpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:2',err
endif
allocate(p%baseprecv(2)%dprcparm(dfpsz),stat=err)
if (err/=0) then
write(0,*)'Precset Memory Failure 2l:3',err
endif
p%baseprecv(2)%iprcparm(:) = 0
p%baseprecv(2)%iprcparm(p_type_) = bja_
p%baseprecv(2)%iprcparm(restr_) = psb_none_
p%baseprecv(2)%iprcparm(prol_) = psb_none_
p%baseprecv(2)%iprcparm(iren_) = 0
p%baseprecv(2)%iprcparm(n_ovr_) = 0
p%baseprecv(2)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(2)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(2)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(2)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(2)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(2)%iprcparm(glb_smth_) = 1
p%baseprecv(2)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(2)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(2)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(2)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(2)%iprcparm(jac_sweeps_) = 1
!!$ allocate(p%baseprecv(ilev_)%iprcparm(ifpsz),stat=err)
!!$ if (err/=0) then
!!$ write(0,*)'Precset Memory Failure 2l:2',err
!!$ endif
!!$ allocate(p%baseprecv(ilev_)%dprcparm(dfpsz),stat=err)
!!$ if (err/=0) then
!!$ write(0,*)'Precset Memory Failure 2l:3',err
!!$ endif
p%baseprecv(ilev_)%iprcparm(:) = 0
p%baseprecv(ilev_)%iprcparm(p_type_) = bja_
p%baseprecv(ilev_)%iprcparm(restr_) = psb_none_
p%baseprecv(ilev_)%iprcparm(prol_) = psb_none_
p%baseprecv(ilev_)%iprcparm(iren_) = 0
p%baseprecv(ilev_)%iprcparm(n_ovr_) = 0
p%baseprecv(ilev_)%iprcparm(ml_type_) = mult_ml_prec_
p%baseprecv(ilev_)%iprcparm(aggr_alg_) = loc_aggr_
p%baseprecv(ilev_)%iprcparm(smth_kind_) = smth_omg_
p%baseprecv(ilev_)%iprcparm(coarse_mat_) = mat_distr_
p%baseprecv(ilev_)%iprcparm(smth_pos_) = post_smooth_
p%baseprecv(ilev_)%iprcparm(glb_smth_) = 1
p%baseprecv(ilev_)%iprcparm(om_choice_) = lib_choice_
p%baseprecv(ilev_)%iprcparm(f_type_) = f_ilu_n_
p%baseprecv(ilev_)%iprcparm(ilu_fill_in_) = 0
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = 4.d0/3.d0
p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = 1
if (present(iv)) then
isz = size(iv)
if (isz >= 1) p%baseprecv(2)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(2)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(2)%iprcparm(smth_kind_) = iv(3)
if (isz >= 4) p%baseprecv(2)%iprcparm(coarse_mat_) = iv(4)
if (isz >= 5) p%baseprecv(2)%iprcparm(smth_pos_) = iv(5)
if (isz >= 6) p%baseprecv(2)%iprcparm(f_type_) = iv(6)
if (isz >= 7) p%baseprecv(2)%iprcparm(jac_sweeps_) = iv(7)
if (isz >= 1) p%baseprecv(ilev_)%iprcparm(ml_type_) = iv(1)
if (isz >= 2) p%baseprecv(ilev_)%iprcparm(aggr_alg_) = iv(2)
if (isz >= 3) p%baseprecv(ilev_)%iprcparm(smth_kind_) = iv(3)
if (isz >= 4) p%baseprecv(ilev_)%iprcparm(coarse_mat_) = iv(4)
if (isz >= 5) p%baseprecv(ilev_)%iprcparm(smth_pos_) = iv(5)
if (isz >= 6) p%baseprecv(ilev_)%iprcparm(f_type_) = iv(6)
if (isz >= 7) p%baseprecv(ilev_)%iprcparm(jac_sweeps_) = iv(7)
end if
if (present(rs)) then
p%baseprecv(2)%iprcparm(om_choice_) = user_choice_
p%baseprecv(2)%dprcparm(smooth_omega_) = rs
p%baseprecv(ilev_)%iprcparm(om_choice_) = user_choice_
p%baseprecv(ilev_)%dprcparm(smooth_omega_) = rs
end if

67
src/prec/psb_zbaseprc_aply.f90
Unescape Escape View File

@ -34,9 +34,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
@ -52,7 +52,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type),intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: prec
complex(kind(0.d0)),intent(inout) :: x(:), y(:)
complex(kind(0.d0)),intent(in) :: beta
complex(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
complex(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -68,13 +68,13 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name, ch_err
interface psb_bjac_aply
subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type), intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: prec
complex(kind(0.d0)),intent(inout) :: x(:), y(:)
complex(kind(0.d0)),intent(in) :: beta
complex(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
complex(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -105,33 +105,35 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
case(noprec_)
n_row=desc_data%matrix_data(psb_n_row_)
if (beta==zzero) then
y(1:n_row) = x(1:n_row)
else if (beta==zone) then
y(1:n_row) = x(1:n_row) + y(1:n_row)
else if (beta==-zone) then
y(1:n_row) = x(1:n_row) - y(1:n_row)
else
y(1:n_row) = x(1:n_row) + beta*y(1:n_row)
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(diagsc_)
n_row=desc_data%matrix_data(psb_n_row_)
if (beta==zzero) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row)
else if (beta==zone) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + y(1:n_row)
else if (beta==-zone) then
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) - y(1:n_row)
if (size(work) >= size(x)) then
ww => work
else
y(1:n_row) = x(1:n_row)*prec%d(1:n_row) + beta*y(1:n_row)
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
end if
n_row=desc_data%matrix_data(psb_n_row_)
ww(1:n_row) = x(1:n_row)*prec%d(1:n_row)
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (size(work) < size(x)) then
deallocate(ww,stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Deallocate')
goto 9999
end if
end if
case(bja_)
call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjac_aply'
@ -142,7 +144,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
if (prec%iprcparm(n_ovr_)==0) then
! shortcut: this fixes performance for RAS(0) == BJA
call psb_bjac_aply(prec,x,beta,y,desc_data,trans,work,info)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjacaply'
@ -214,7 +216,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
end if
endif
call psb_bjac_aply(prec,tx,zzero,ty,prec%desc_data,trans,aux,info)
call psb_bjac_aply(zone,prec,tx,zzero,ty,prec%desc_data,trans,aux,info)
if(info.ne.0) then
info=4010
ch_err='psb_bjac_aply'
@ -250,18 +252,7 @@ subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
& prec%iprcparm(prol_)
end select
if (beta == zzero) then
y(1:desc_data%matrix_data(psb_n_row_)) = ty(1:desc_data%matrix_data(psb_n_row_))
else if (beta == zone) then
y(1:desc_data%matrix_data(psb_n_row_)) = y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
else if (beta == -zone) then
y(1:desc_data%matrix_data(psb_n_row_)) = -y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
else
y(1:desc_data%matrix_data(psb_n_row_)) = beta*y(1:desc_data%matrix_data(psb_n_row_)) + &
& ty(1:desc_data%matrix_data(psb_n_row_))
end if
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
if ((6*isz) <= size(work)) then

15
src/prec/psb_zbaseprc_bld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd)
Implicit None
type(psb_zspmat_type), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_zbaseprc_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -88,7 +88,7 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd)
use psb_const_mod
implicit none
type(psb_zspmat_type), intent(in) :: a
type(psb_zspmat_type), intent(inout) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_zbaseprc_type), intent(inout) :: p
integer, intent(out) :: info
@ -169,7 +169,13 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd)
select case(p%iprcparm(p_type_))
case (noprec_)
! Do nothing.
call psb_cdcpy(desc_a,p%desc_data,info)
if(info /= 0) then
info=4010
ch_err='psb_cdcpy'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (diagsc_)
@ -256,7 +262,8 @@ subroutine psb_zbaseprc_bld(a,desc_a,p,info,upd)
end select
p%base_a => a
p%base_desc => desc_a
call psb_erractionrestore(err_act)
return

41
src/prec/psb_zbjac_aply.f90
Unescape Escape View File

@ -34,9 +34,9 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a Block Jacobi preconditioner stored in prec
! Note that desc_data may or may not be the same as prec%desc_data,
! but since both are INTENT(IN) this should be legal.
@ -54,7 +54,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type), intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: prec
complex(kind(0.d0)),intent(inout) :: x(:), y(:)
complex(kind(0.d0)),intent(in) :: beta
complex(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
complex(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -125,7 +125,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
& trans='N',unit=diagl,choice=psb_none_,work=aux)
if(info /=0) goto 9999
ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row)
call psb_spsm(zone,prec%av(u_pr_),ww,beta,y,desc_data,info,&
call psb_spsm(alpha,prec%av(u_pr_),ww,beta,y,desc_data,info,&
& trans='N',unit=diagu,choice=psb_none_, work=aux)
if(info /=0) goto 9999
@ -134,7 +134,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
& trans=trans,unit=diagu,choice=psb_none_, work=aux)
if(info /=0) goto 9999
ww(1:n_row) = ww(1:n_row)*prec%d(1:n_row)
call psb_spsm(zone,prec%av(l_pr_),ww,beta,y,desc_data,info,&
call psb_spsm(alpha,prec%av(l_pr_),ww,beta,y,desc_data,info,&
& trans=trans,unit=diagl,choice=psb_none_,work=aux)
if(info /=0) goto 9999
@ -152,16 +152,8 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
end select
if(info /=0) goto 9999
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (beta == zzero) then
y(1:n_row) = ww(1:n_row)
else if (beta==zone) then
y(1:n_row) = ww(1:n_row) + y(1:n_row)
else if (beta==-zone) then
y(1:n_row) = ww(1:n_row) - y(1:n_row)
else
y(1:n_row) = ww(1:n_row) + beta*y(1:n_row)
endif
case (f_umf_)
@ -174,15 +166,7 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
if (beta == zzero) then
y(1:n_row) = ww(1:n_row)
else if (beta==zone) then
y(1:n_row) = ww(1:n_row) + y(1:n_row)
else if (beta==-zone) then
y(1:n_row) = ww(1:n_row) - y(1:n_row)
else
y(1:n_row) = ww(1:n_row) + beta*y(1:n_row)
endif
call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
case default
write(0,*) 'Unknown factorization type in bjac_aply',prec%iprcparm(f_type_)
@ -253,15 +237,8 @@ subroutine psb_zbjac_aply(prec,x,beta,y,desc_data,trans,work,info)
end select
if (beta == zzero) then
y(1:n_row) = tx(1:n_row)
else if (beta==zone) then
y(1:n_row) = tx(1:n_row) + y(1:n_row)
else if (beta==-zone) then
y(1:n_row) = tx(1:n_row) - y(1:n_row)
else
y(1:n_row) = tx(1:n_row) + beta*y(1:n_row)
endif
call psb_geaxpby(alpha,tx,beta,y,desc_data,info)
deallocate(tx,ty)

21
src/prec/psb_zbldaggrmat.f90
Unescape Escape View File

@ -620,7 +620,12 @@ contains
! Doing it this way means to consider diag(Ai)
!
!
call psb_symbmm(am3,am4,am1)
call psb_symbmm(am3,am4,am1,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
goto 9999
end if
call psb_numbmm(am3,am4,am1)
if (debug) write(0,*) me,'Done NUMBMM 1'
@ -667,7 +672,12 @@ contains
if (test_dump) &
& call psb_csprt(60+me,am1,head='% (I-wDA)Pt',ivr=desc_a%loc_to_glob)
call psb_symbmm(a,am1,am3)
call psb_symbmm(a,am1,am3,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
goto 9999
end if
call psb_numbmm(a,am1,am3)
if (debug) write(0,*) me,'Done NUMBMM 2'
@ -724,7 +734,12 @@ contains
endif
if (debug) write(0,*) me,'starting symbmm 3'
call psb_symbmm(am2,am3,b)
call psb_symbmm(am2,am3,b,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='symbmm')
goto 9999
end if
if (debug) write(0,*) me,'starting numbmm 3'
call psb_numbmm(am2,am3,b)
if (debug) write(0,*) me,'Done NUMBMM 3'

687
src/prec/psb_zmlprc_aply.f90
Unescape Escape View File

@ -34,10 +34,53 @@
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + K^-1 X
! where K is a multilevel (actually 2-level) preconditioner stored in prec
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a multilevel preconditioner stored in baseprecv
!
! cfr.: Smith, Biorstad & Gropp
! Domain Decomposition
! Cambridge Univ. Press
!
! To each level I there corresponds a matrix A(I) and a preconditioner K(I)
!
! A notational difference: in the DD reference above the preconditioner for
! a given level K(I) is written out as a sum over the subdomains
!
! SUM_k(R_k^T A_k R_k)
!
! whereas in this code the sum is implicit in the parallelization,
! i.e. each process takes care of one subdomain, and for each level we have
! as many subdomains as there are processes (except for the coarsest level where
! we might have a replicated index space). Thus the sum apparently disappears
! from our code, but only apparently, because it is implicit in the call
! to psb_baseprc_aply.
!
! A bit of description of the baseprecv(:) data structure:
! 1. Number of levels = NLEV = size(baseprecv(:))
! 2. baseprecv(ilev)%av(:) sparse matrices needed for the current level.
! Includes:
! 2.1.: baseprecv(ilev)%av(l_pr_) L factor of ILU preconditioners
! 2.2.: baseprecv(ilev)%av(u_pr_) U factor of ILU preconditioners
! 2.3.: baseprecv(ilev)%av(ap_nd_) Off-diagonal part of A for Jacobi sweeps
! 2.4.: baseprecv(ilev)%av(ac_) Aggregated matrix of level ILEV
! 2.5.: baseprecv(ilev)%av(sm_pr_t_) Smoother prolongator transpose; maps vectors
! (ilev-1) ---> (ilev)
! 2.6.: baseprecv(ilev)%av(sm_pr_) Smoother prolongator; maps vectors
! (ilev) ---> (ilev-1)
! Shouldn't we keep just one of them and handle transpose in the sparse BLAS? maybe
!
! 3. baseprecv(ilev)%desc_data comm descriptor for level ILEV
! 4. baseprecv(ilev)%base_a Pointer (really a pointer!) to the base matrix
! of the current level, i.e.: if ILEV=1 then A
! else the aggregated matrix av(ac_); so we have
! a unified treatment of residuals. Need this to
! avoid passing explicitly matrix A to the
! outer prec. routine
! 5. baseprecv(ilev)%mlia The aggregation map from (ilev-1)-->(ilev)
! if no smoother, it is used instead of sm_pr_
! 6. baseprecv(ilev)%nlaggr Number of aggregates on the various procs.
!
use psb_serial_mod
@ -52,7 +95,7 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
type(psb_desc_type),intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: baseprecv(:)
complex(kind(1.d0)),intent(in) :: beta
complex(kind(1.d0)),intent(in) :: alpha,beta
complex(kind(1.d0)),intent(inout) :: x(:), y(:)
character :: trans
complex(kind(1.d0)),target :: work(:)
@ -68,25 +111,32 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
real(kind(1.d0)) :: omega
real(kind(1.d0)) :: t1, t2, t3, t4, t5, t6, t7, mpi_wtime
logical, parameter :: debug=.false., debugprt=.false.
integer :: ismth
integer :: ismth, nlev, ilev
external mpi_wtime
character(len=20) :: name, ch_err
type psb_mlprec_wrk_type
complex(kind(1.d0)), pointer :: tx(:)=>null(),ty(:)=>null(),&
& x2l(:)=>null(),y2l(:)=>null(),&
& b2l(:)=>null(),tty(:)=>null()
end type psb_mlprec_wrk_type
type(psb_mlprec_wrk_type), pointer :: mlprec_wrk(:)
interface psb_baseprc_aply
subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: prec
complex(kind(1.d0)),intent(inout) :: x(:), y(:)
complex(kind(1.d0)),intent(in) :: beta
complex(kind(1.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
complex(kind(1.d0)),target :: work(:)
integer, intent(out) :: info
end subroutine psb_zbaseprc_aply
end interface
name='psb_zmlprc_aply'
name='psb_mlprc_aply'
info = 0
call psb_erractionsave(err_act)
@ -94,426 +144,625 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
omega=baseprecv(2)%dprcparm(smooth_omega_)
ismth=baseprecv(2)%iprcparm(smth_kind_)
nlev = size(baseprecv)
allocate(mlprec_wrk(nlev),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
select case(baseprecv(2)%iprcparm(ml_type_))
case(no_ml_)
! Should not really get here.
write(0,*) 'Smooth preconditioner with no multilevel in MLPRC_APLY????'
call psb_errpush(4010,name,a_err='no_ml_ in mlprc_aply?')
goto 9999
case(add_ml_prec_)
!
! Additive multilevel
! Additive is very simple.
! 1. X(1) = Xext
! 2. DO ILEV=2,NLEV
! X(ILEV) = AV(PR_SM_T_)*X(ILEV-1)
! 3. Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! 4. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = AV(PR_SM_)*Y(ILEV+1)
! 5. Yext = beta*Yext + Y(1)
!
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
call psb_baseprc_aply(baseprecv(1),x,beta,y,desc_data,trans,work,info)
if(info /=0) goto 9999
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),stat=info)
call psb_baseprc_aply(alpha,baseprecv(1),x,beta,y,&
& baseprecv(1)%base_desc,trans,work,info)
if(info /=0) goto 9999
allocate(mlprec_wrk(1)%x2l(size(x)),mlprec_wrk(1)%y2l(size(y)))
mlprec_wrk(1)%x2l(:) = x(:)
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
allocate(mlprec_wrk(ilev)%x2l(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%tx(max(n_row,n_col)),&
& mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = zzero
w2l(:) = zzero
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(1:n_row) = mlprec_wrk(ilev-1)%x2l(1:n_row)
mlprec_wrk(ilev)%tx(n_row+1:max(n_row,n_col)) = zzero
mlprec_wrk(ilev)%ty(:) = zzero
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
!
! Smoothed aggregation
!
allocate(tx(max(n_row,n_col)),ty(max(n_row,n_col)),&
& tz(max(n_row,n_col)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
tx(1:desc_data%matrix_data(psb_n_row_)) = x(1:desc_data%matrix_data(psb_n_row_))
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
ty(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
tz(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = zzero
end if
call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info)
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcut
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + x(i)
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%x2l(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_baseprc_aply(zone,baseprecv(ilev),&
& mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev)%desc_data, 'N',work,info)
enddo
do ilev =nlev,2,-1
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
if (ismth /= no_smth_) then
call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info)
if(info /=0) goto 9999
!
! Finally add back into Y.
!
call psb_geaxpby(zone,ty,zone,y,desc_data,info)
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_),mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,info)
if(info /=0) goto 9999
deallocate(tx,ty,tz)
else
do i=1, desc_data%matrix_data(psb_n_row_)
y(i) = y(i) + t2l(baseprecv(2)%mlia(i))
do i=1, n_row
mlprec_wrk(ilev-1)%y2l(i) = mlprec_wrk(ilev-1)%y2l(i) + &
& mlprec_wrk(ilev)%y2l(baseprecv(ilev)%mlia(i))
enddo
end if
end do
if (debugprt) write(0,*)' Y2: ',Y(:)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,zone,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
deallocate(t2l,w2l)
case(mult_ml_prec_)
!
! Multiplicative multilevel
! Pre/post smoothing versions.
!
select case(baseprecv(2)%iprcparm(smth_pos_))
case(post_smooth_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = zzero
w2l(:) = zzero
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==zzero
! Post smoothing.
! 1. X(1) = Xext
! 2. DO ILEV=2, NLEV :: X(ILEV) = AV(PR_SM_T_,ILEV)*X(ILEV-1)
! 3. Y(NLEV) = (K(NLEV)**(-1))*X(NLEV)
! 4. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV))
!
! 5. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
! Also: post smoothing is not spelled out in detail in DD.
!
!
if (debug) write(0,*)' mult_ml_apply omega ',omega
if (debugprt) write(0,*)' mult_ml_apply X: ',X(:)
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%tx(nr2l), stat=info)
mlprec_wrk(1)%x2l(:) = zzero
mlprec_wrk(1)%y2l(:) = zzero
mlprec_wrk(1)%tx(:) = zzero
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,&
& baseprecv(1)%base_desc,info)
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
do ilev=2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth = baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
if (info /= 0) then
if (debug) write(0,*)' From axpby1 ',size(x),size(tx),n_row,n_col,nr2l,nrg
call psb_errpush(4010,name,a_err='axpby post_smooth 1')
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= no_smth_) then
!
! Smoothed aggregation
!
allocate(tz(max(n_row,n_col)),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
mlprec_wrk(ilev-1)%x2l(n_row+1:max(n_row,n_col)) = zzero
end if
call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info)
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%x2l, &
& zzero,mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcut
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%x2l(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) Then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
enddo
call psb_baseprc_aply(zone,baseprecv(nlev),mlprec_wrk(nlev)%x2l, &
& zzero, mlprec_wrk(nlev)%y2l,baseprecv(nlev)%desc_data,'N',work,info)
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,&
&'N',work,info)
if(info /=0) goto 9999
do ilev=nlev-1, 1, -1
ismth = baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info)
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zzero,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
!
! Finally add back into Y.
!
deallocate(tz)
else
ty(:) = zzero
do i=1, desc_data%matrix_data(psb_n_row_)
ty(i) = ty(i) + t2l(baseprecv(2)%mlia(i))
n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_)
mlprec_wrk(ilev)%y2l(:) = zzero
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
deallocate(t2l,w2l)
call psb_spmm(-zone,baseprecv(2)%aorig,ty,zone,tx,desc_data,info,work=work)
if(info /=0) goto 9999
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
call psb_baseprc_aply(baseprecv(1),tx,zone,ty,desc_data,trans,&
& work,info)
if(info /=0) goto 9999
call psb_geaxpby(zone,ty,beta,y,desc_data,info)
call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
deallocate(tx,ty)
enddo
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
case(pre_smooth_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info)
!
! Pre smoothing.
! 1. X(1) = Xext
! 2. Y(1) = (K(1)**(-1))*X(1)
! 3. TX(1) = X(1) - A(1)*Y(1)
! 4. DO ILEV=2, NLEV
! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1)
! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV))
! 5. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! 6. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
!
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%tx(nr2l), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = zzero
w2l(:) = zzero
mlprec_wrk(1)%y2l(:) = zzero
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==zero
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if(info /=0) goto 9999
call psb_baseprc_aply(baseprecv(1),x,zzero,tty,desc_data,&
mlprec_wrk(1)%x2l(:) = x
call psb_baseprc_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,&
& zzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
if(info /=0) goto 9999
call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work)
mlprec_wrk(1)%tx = mlprec_wrk(1)%x2l
call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%tx,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
if (ismth /= no_smth_) then
allocate(tz(max(n_row,n_col)),stat=info)
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth = baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%tx(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,desc_data,info,work=work)
if (ismth /= no_smth_) then
!
!Smoothed Aggregation
!
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
mlprec_wrk(ilev-1)%tx(n_row+1:max(n_row,n_col)) = zzero
end if
call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info)
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%tx,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcuts
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
mlprec_wrk(ilev)%x2l = zzero
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%tx(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if(info /=0) goto 9999
if(ilev < nlev) then
mlprec_wrk(ilev)%tx = mlprec_wrk(ilev)%x2l
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
endif
enddo
do ilev = nlev-1, 1, -1
ismth=baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info)
if(info /=0) goto 9999
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,info)
call psb_geaxpby(zone,ty,zone,tty,desc_data,info)
if(info /=0) goto 9999
deallocate(tz)
else
do i=1, desc_data%matrix_data(psb_n_row_)
tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i))
n_row = baseprecv(ilev+1)%base_desc%matrix_data(psb_n_row_)
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
call psb_geaxpby(zone,tty,beta,y,desc_data,info)
enddo
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
if(info /=0) goto 9999
deallocate(t2l,w2l,tx,ty,tty)
case(smooth_both_)
t1 = mpi_wtime()
n_row = desc_data%matrix_data(psb_n_row_)
n_col = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(2)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(2)%desc_data%matrix_data(psb_n_row_)
allocate(t2l(nr2l),w2l(nr2l),tx(n_col),ty(n_col),tty(n_col),stat=info)
!
! Symmetrized smoothing.
! 1. X(1) = Xext
! 2. Y(1) = (K(1)**(-1))*X(1)
! 3. TX(1) = X(1) - A(1)*Y(1)
! 4. DO ILEV=2, NLEV
! X(ILEV) = AV(PR_SM_T_,ILEV)*TX(ILEV-1)
! Y(ILEV) = (K(ILEV)**(-1))*X(ILEV)
! TX(ILEV) = (X(ILEV)-A(ILEV)*Y(ILEV))
! 5. DO ILEV=NLEV-1,1,-1
! Y(ILEV) = Y(ILEV) + AV(PR_SM_,ILEV+1)*Y(ILEV+1)
! Y(ILEV) = Y(ILEV) + (K(ILEV)**(-1))*(X(ILEV)-A(ILEV)*Y(ILEV))
! 6. Yext = beta*Yext + Y(1)
!
! Note: level numbering reversed wrt ref. DD, i.e.
! 1..NLEV <=> (j) <-> 0
!
!
n_col = desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(1)%desc_data%matrix_data(psb_n_col_)
allocate(mlprec_wrk(1)%x2l(nr2l),mlprec_wrk(1)%y2l(nr2l), &
& mlprec_wrk(1)%ty(nr2l), mlprec_wrk(1)%tx(nr2l), stat=info)
mlprec_wrk(1)%x2l(:) = zzero
mlprec_wrk(1)%y2l(:) = zzero
mlprec_wrk(1)%tx(:) = zzero
mlprec_wrk(1)%ty(:) = zzero
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
t2l(:) = zzero
w2l(:) = zzero
tx(:) = zzero
ty(:) = zzero
tty(:) = zzero
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%tx,&
& baseprecv(1)%base_desc,info)
!
! Need temp copies to handle Y<- betaY + K^-1 X
! One of the temp copies is not strictly needed when beta==zero
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if(info /=0) goto 9999
call psb_baseprc_aply(zone,baseprecv(1),mlprec_wrk(1)%x2l,&
& zzero,mlprec_wrk(1)%y2l,&
& baseprecv(1)%base_desc,&
& trans,work,info)
call psb_baseprc_aply(baseprecv(1),tx,zzero,tty,desc_data,trans,work,info)
if(info /=0) goto 9999
call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work)
mlprec_wrk(1)%ty = mlprec_wrk(1)%x2l
call psb_spmm(-zone,baseprecv(1)%base_a,mlprec_wrk(1)%y2l,&
& zone,mlprec_wrk(1)%ty,baseprecv(1)%base_desc,info,work=work)
if(info /=0) goto 9999
do ilev = 2, nlev
n_row = baseprecv(ilev-1)%base_desc%matrix_data(psb_n_row_)
n_col = baseprecv(ilev-1)%desc_data%matrix_data(psb_n_col_)
nr2l = baseprecv(ilev)%desc_data%matrix_data(psb_n_col_)
nrg = baseprecv(ilev)%desc_data%matrix_data(psb_n_row_)
ismth=baseprecv(ilev)%iprcparm(smth_kind_)
allocate(mlprec_wrk(ilev)%ty(nr2l),mlprec_wrk(ilev)%y2l(nr2l),&
& mlprec_wrk(ilev)%x2l(nr2l), stat=info)
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
mlprec_wrk(ilev)%ty(:) = zzero
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
if (ismth /= no_smth_) then
if (baseprecv(2)%iprcparm(glb_smth_) >0) then
call psb_halo(tx,baseprecv(1)%desc_data,info,work=work)
!
!Smoothed Aggregation
!
if (baseprecv(ilev)%iprcparm(glb_smth_) >0) then
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if(info /=0) goto 9999
else
tx(desc_data%matrix_data(psb_n_row_)+1:max(n_row,n_col)) = zzero
mlprec_wrk(ilev-1)%ty(n_row+1:max(n_row,n_col)) = zzero
end if
call psb_csmm(zone,baseprecv(2)%av(sm_pr_t_),tx,zzero,t2l,info)
call psb_csmm(zone,baseprecv(ilev)%av(sm_pr_t_),mlprec_wrk(ilev-1)%ty,zzero,&
& mlprec_wrk(ilev)%x2l,info)
if(info /=0) goto 9999
else
!
! Raw aggregation, may take shortcuts
!
do i=1,desc_data%matrix_data(psb_n_row_)
t2l(baseprecv(2)%mlia(i)) = t2l(baseprecv(2)%mlia(i)) + tx(i)
mlprec_wrk(ilev)%x2l = zzero
do i=1,n_row
mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) = &
& mlprec_wrk(ilev)%x2l(baseprecv(ilev)%mlia(i)) + &
& mlprec_wrk(ilev-1)%ty(i)
end do
end if
if (baseprecv(2)%iprcparm(coarse_mat_)==mat_repl_) Then
call psb_sum(ictxt,t2l(1:nrg))
else if (baseprecv(2)%iprcparm(coarse_mat_) /= mat_distr_) Then
if (baseprecv(ilev)%iprcparm(coarse_mat_)==mat_repl_) then
call psb_sum(ictxt,mlprec_wrk(ilev)%x2l(1:nrg))
else if (baseprecv(ilev)%iprcparm(coarse_mat_) /= mat_distr_) then
write(0,*) 'Unknown value for baseprecv(2)%iprcparm(coarse_mat_) ',&
& baseprecv(2)%iprcparm(coarse_mat_)
& baseprecv(ilev)%iprcparm(coarse_mat_)
endif
t6 = mpi_wtime()
w2l=t2l
call psb_baseprc_aply(baseprecv(2),w2l,zzero,t2l,baseprecv(2)%desc_data,&
& 'N',work,info)
call psb_geaxpby(zone,mlprec_wrk(ilev)%x2l,zzero,mlprec_wrk(ilev)%tx,&
& baseprecv(ilev)%base_desc,info)
if(info /=0) goto 9999
if (ismth /= no_smth_) then
call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%x2l,&
& zzero,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%desc_data, 'N',work,info)
if (ismth == smth_omg_) &
& call psb_halo(t2l,baseprecv(2)%desc_data,info,work=work)
call psb_csmm(zone,baseprecv(2)%av(sm_pr_),t2l,zzero,ty,info)
if(info /=0) goto 9999
call psb_geaxpby(zone,ty,zone,tty,desc_data,info)
if(ilev < nlev) then
mlprec_wrk(ilev)%ty = mlprec_wrk(ilev)%x2l
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%ty,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
endif
else
enddo
do i=1, desc_data%matrix_data(psb_n_row_)
tty(i) = tty(i) + t2l(baseprecv(2)%mlia(i))
do ilev=nlev-1, 1, -1
ismth=baseprecv(ilev+1)%iprcparm(smth_kind_)
if (ismth /= no_smth_) then
if (ismth == smth_omg_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%desc_data,&
& info,work=work)
call psb_csmm(zone,baseprecv(ilev+1)%av(sm_pr_),mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,info)
if(info /=0) goto 9999
else
n_row = baseprecv(ilev)%base_desc%matrix_data(psb_n_row_)
do i=1, n_row
mlprec_wrk(ilev)%y2l(i) = mlprec_wrk(ilev)%y2l(i) + &
& mlprec_wrk(ilev+1)%y2l(baseprecv(ilev+1)%mlia(i))
enddo
end if
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_spmm(-zone,baseprecv(ilev)%base_a,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev)%tx,baseprecv(ilev)%base_desc,info,work=work)
if(info /=0) goto 9999
call psb_spmm(-zone,baseprecv(2)%aorig,tty,zone,tx,desc_data,info,work=work)
call psb_baseprc_aply(zone,baseprecv(ilev),mlprec_wrk(ilev)%tx,&
& zone,mlprec_wrk(ilev)%y2l,baseprecv(ilev)%base_desc, trans, work,info)
if(info /=0) goto 9999
call psb_baseprc_aply(baseprecv(1),tx,zone,tty,desc_data,'N',work,info)
enddo
call psb_geaxpby(zone,tty,beta,y,desc_data,info)
call psb_geaxpby(alpha,mlprec_wrk(1)%y2l,beta,y,&
& baseprecv(1)%base_desc,info)
deallocate(t2l,w2l,tx,ty,tty)
if(info /=0) goto 9999
case default
write(0,*) 'Unknown value for ml_smooth_pos',baseprecv(2)%iprcparm(smth_pos_)
call psb_errpush(4013,name,a_err='wrong smooth_pos',&
& i_Err=(/baseprecv(2)%iprcparm(smth_pos_),0,0,0,0/))
goto 9999
end select
case default
write(0,*) me, 'Wrong mltype into PRC_APLY ',&
& baseprecv(2)%iprcparm(ml_type_)
call psb_errpush(4013,name,a_err='wrong mltype',&
& i_Err=(/baseprecv(2)%iprcparm(ml_type_),0,0,0,0/))
goto 9999
end select
call mlprec_wrk_free(mlprec_wrk)
deallocate(mlprec_wrk)
call psb_erractionrestore(err_act)
return
@ -526,4 +775,24 @@ subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
end if
return
contains
subroutine mlprec_wrk_free(wrk)
type(psb_mlprec_wrk_type) :: wrk(:)
! This will not be needed when we have allocatables, as
! it is sufficient to deallocate the container, and
! the compiler is supposed to recursively deallocate the
! various components.
integer i
do i=1, size(wrk)
if (associated(wrk(i)%tx)) deallocate(wrk(i)%tx)
if (associated(wrk(i)%ty)) deallocate(wrk(i)%ty)
if (associated(wrk(i)%x2l)) deallocate(wrk(i)%x2l)
if (associated(wrk(i)%y2l)) deallocate(wrk(i)%y2l)
if (associated(wrk(i)%b2l)) deallocate(wrk(i)%b2l)
if (associated(wrk(i)%tty)) deallocate(wrk(i)%tty)
end do
end subroutine mlprec_wrk_free
end subroutine psb_zmlprc_aply

6
src/prec/psb_zmlprc_bld.f90
Unescape Escape View File

@ -135,8 +135,6 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
p%aorig => a
nullify(p%d)
@ -180,8 +178,8 @@ subroutine psb_zmlprc_bld(a,desc_a,p,info)
! We have used a separate ac because:
! 1. We want to reuse the same routines psb_ilu_bld etc.
! 2. We do NOT want to pass an argument twice to them
! p%av(ac_) and p
! Hence a separate AC and a TRANSFER function.
! p%av(ac_) and p, as this would violate the Fortran standard
! Hence a separate AC and a TRANSFER function at the end.
!
call psb_sp_transfer(ac,p%av(ac_),info)

12
src/prec/psb_zprc_aply.f90
Unescape Escape View File

@ -61,13 +61,13 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
character(len=20) :: name
interface psb_baseprc_aply
subroutine psb_zbaseprc_aply(prec,x,beta,y,desc_data,trans,work,info)
subroutine psb_zbaseprc_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: prec
complex(kind(0.d0)),intent(inout) :: x(:), y(:)
complex(kind(0.d0)),intent(in) :: beta
complex(kind(0.d0)),intent(in) :: alpha,beta
character(len=1) :: trans
complex(kind(0.d0)),target :: work(:)
integer, intent(out) :: info
@ -75,12 +75,12 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
end interface
interface psb_mlprc_aply
subroutine psb_zmlprc_aply(baseprecv,x,beta,y,desc_data,trans,work,info)
subroutine psb_zmlprc_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
use psb_descriptor_type
use psb_prec_type
type(psb_desc_type),intent(in) :: desc_data
type(psb_zbaseprc_type), intent(in) :: baseprecv(:)
complex(kind(0.d0)),intent(in) :: beta
complex(kind(0.d0)),intent(in) :: alpha,beta
complex(kind(0.d0)),intent(inout) :: x(:), y(:)
character :: trans
complex(kind(0.d0)),target :: work(:)
@ -117,14 +117,14 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
end if
if (size(prec%baseprecv) >1) then
if (debug) write(0,*) 'Into mlprc_aply',size(x),size(y)
call psb_mlprc_aply(prec%baseprecv,x,zzero,y,desc_data,trans_,work_,info)
call psb_mlprc_aply(zone,prec%baseprecv,x,zzero,y,desc_data,trans_,work_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_zmlprc_aply')
goto 9999
end if
else if (size(prec%baseprecv) == 1) then
call psb_baseprc_aply(prec%baseprecv(1),x,zzero,y,desc_data,trans_, work_,info)
call psb_baseprc_aply(zone,prec%baseprecv(1),x,zzero,y,desc_data,trans_, work_,info)
else
write(0,*) 'Inconsistent preconditioner: size of baseprecv???'
endif

32
src/prec/psb_zprecbld.f90
Unescape Escape View File

@ -49,7 +49,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
Implicit None
type(psb_zspmat_type), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_zprec_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -60,7 +60,7 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
use psb_descriptor_type
use psb_prec_type
type(psb_zspmat_type), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_desc_type), intent(in), target :: desc_a
type(psb_zbaseprc_type),intent(inout) :: p
integer, intent(out) :: info
character, intent(in), optional :: upd
@ -149,7 +149,9 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
endif
if (size(p%baseprecv) > 1) then
call init_baseprc_av(p%baseprecv(2),info)
do i=2, size(p%baseprecv)
call init_baseprc_av(p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='allocate'
@ -157,31 +159,17 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
goto 9999
endif
call psb_mlprc_bld(a,desc_a,p%baseprecv(2),info)
call psb_mlprc_bld(p%baseprecv(i-1)%base_a,p%baseprecv(i-1)%base_desc,&
& p%baseprecv(i),info)
if (info /= 0) then
info=4010
ch_err='psb_mlprc_bld'
call psb_errpush(info,name,a_err=ch_err)
call psb_errpush(info,name)
goto 9999
endif
!
! 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
if (debug) then
write(0,*) 'Return from ',i-1,' call to mlprcbld ',info
endif
call psb_mlprc_bld(p%baseprecv(i-1)%av(ac_),p%baseprecv(i-1)%desc_data,&
& p%baseprecv(i),info)
end do
endif

2
src/serial/dp/check_dim.f
Unescape Escape View File

@ -51,7 +51,7 @@ C Local scalars
NNZ = NZ
LIMIT = INT(DIM_BLOCK*PSB_PERCENT_)
c$$$ LIMIT = INT(DIM_BLOCK*PSB_PERCENT_)
DO BLOCK = 1, NG
DIM_BLOCK = IA(1,BLOCK+1)-IA(1,BLOCK)

17
src/serial/f77/smmp.f
Unescape Escape View File

@ -24,10 +24,21 @@ c
* index(*)
integer, pointer :: ic(:),jc(:)
integer :: nze, info
integer, save :: iunit=11
c
c symbolic matrix multiply c=a*b
c
c$$$ write(0,*) 'SYMBMM: ',n,m,l,ib(m+1)-1,jb(ib(m+1)-1)
c$$$ open(iunit)
c$$$ write(iunit,*) 'SYMBMM: ',n,m,l
c$$$ write(iunit,*) 'SYMBMM: A:'
c$$$ do i=1,n
c$$$ write(iunit,*) 'Row:',i,' : ',ja(ia(i):ia(i+1)-1)
c$$$ enddo
c$$$ write(iunit,*) 'SYMBMM: B:'
c$$$ do i=1,m
c$$$ write(iunit,*) 'Row:',i,' : ',jb(ib(i):ib(i+1)-1)
c$$$ enddo
if (size(ic) < n+1) then
call psb_realloc(n+1,ic,info)
endif
@ -95,6 +106,7 @@ c
endif
call psb_realloc(nze,jc,info)
end if
do 40 j= ic(i),ic(i+1)-1
if (diagc.eq.1 .and. istart.eq.i) then
istart = index(istart)
@ -105,8 +117,11 @@ c
index(jc(j))=0
40 continue
call isr(length,jc(ic(i)))
c$$$ write(iunit,*) length,' : ',jc(ic(i):ic(i)+length-1)
index(i) = 0
50 continue
c$$$ close(iunit)
c$$$ iunit = iunit + 1
c$$$ write(0,*) 'SYMBMM: on exit',ic(n+1)-1,jc(ic(n+1)-1)
return
end

2
src/serial/jad/djdnrmi.f
Unescape Escape View File

@ -49,7 +49,7 @@ C .. External routines ..
PNG = IA(1)
PIA = IA(2)
PJA = IA(3)
djdnrmi = 0.d0
IF (DESCRA(1:1).EQ.'G') THEN
DJDNRMI = DJADNR(TRANS,M,N,IA(PNG),
+ A,JA,IA(PJA),IA(PIA),

4
src/serial/psb_cest.f90
Unescape Escape View File

@ -66,7 +66,7 @@ subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, iup, info)
lar = nnz
else
info = 136
call psb_errpush(info,name,a_err=afmt)
call psb_errpush(info,name,a_err=toupper(afmt))
goto 9999
endif
@ -86,7 +86,7 @@ subroutine psb_cest(afmt, m,n,nnz, lia1, lia2, lar, iup, info)
lar = nnz
else
info = 136
call psb_errpush(info,name,a_err=afmt)
call psb_errpush(info,name,a_err=toupper(afmt))
goto 9999
endif

4
src/serial/psb_dcsprt.f90
Unescape Escape View File

@ -70,11 +70,12 @@ subroutine psb_dcsprt(iout,a,iv,eirs,eics,head,ivr,ivc)
else
ics = 0
endif
open(iout)
if (present(head)) then
write(iout,'(a)') '%%MatrixMarket matrix coordinate real general'
write(iout,'(a,a)') '% ',head
write(iout,'(a)') '%'
write(iout,'(a,a)') '% ',toupper(a%fida)
endif
select case(toupper(a%fida))
@ -181,4 +182,5 @@ subroutine psb_dcsprt(iout,a,iv,eirs,eics,head,ivr,ivc)
case default
write(0,*) 'Feeling lazy today, format not implemented: "',a%fida,'"'
end select
close(iout)
end subroutine psb_dcsprt

36
src/serial/psb_dsymbmm.f90
Unescape Escape View File

@ -38,8 +38,9 @@
! rewritten in Fortran 95 making use of our sparse matrix facilities.
!
subroutine psb_dsymbmm(a,b,c)
subroutine psb_dsymbmm(a,b,c,info)
use psb_spmat_type
use psb_string_mod
implicit none
type(psb_dspmat_type) :: a,b,c
@ -54,7 +55,6 @@ subroutine psb_dsymbmm(a,b,c)
integer, pointer :: ic(:),jc(:)
end subroutine symbmm
end interface
interface psb_sp_getrow
subroutine psb_dspgetrow(irw,a,nz,ia,ja,val,info,iren,lrw)
use psb_spmat_type
@ -69,6 +69,28 @@ subroutine psb_dsymbmm(a,b,c)
end subroutine psb_dspgetrow
end interface
character(len=20) :: name, ch_err
integer :: err_act
name='psb_symbmm'
call psb_erractionsave(err_act)
select case(toupper(a%fida(1:3)))
case ('CSR')
case default
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
select case(toupper(b%fida(1:3)))
case ('CSR')
case default
info=136
ch_err=b%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
endif
@ -78,7 +100,6 @@ subroutine psb_dsymbmm(a,b,c)
endif
nze = max(a%m+1,2*a%m)
call psb_sp_reall(c,nze,info)
!!$ write(0,*) 'SYMBMM90 ',size(c%pl),size(c%pr)
!
! Note: we need to test whether there is a performance impact
! in not using the original Douglas & Bank code.
@ -97,6 +118,15 @@ subroutine psb_dsymbmm(a,b,c)
c%fida='CSR'
c%descra='GUN'
deallocate(itemp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error()
return
end if
return
contains

33
src/serial/psb_zsymbmm.f90
Unescape Escape View File

@ -38,8 +38,9 @@
! rewritten in Fortran 95 making use of our sparse matrix facilities.
!
subroutine psb_zsymbmm(a,b,c)
subroutine psb_zsymbmm(a,b,c,info)
use psb_spmat_type
use psb_string_mod
implicit none
type(psb_zspmat_type) :: a,b,c
@ -68,6 +69,27 @@ subroutine psb_zsymbmm(a,b,c)
integer, intent(out) :: info
end subroutine psb_zspgetrow
end interface
character(len=20) :: name, ch_err
integer :: err_act
name='psb_symbmm'
call psb_erractionsave(err_act)
select case(toupper(a%fida(1:3)))
case ('CSR')
case default
info=135
ch_err=a%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
select case(toupper(b%fida(1:3)))
case ('CSR')
case default
info=136
ch_err=b%fida(1:3)
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
if (b%m /= a%k) then
write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
@ -97,6 +119,15 @@ subroutine psb_zsymbmm(a,b,c)
c%fida='CSR'
c%descra='GUN'
deallocate(itemp)
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.act_abort) then
call psb_error()
return
end if
return
contains

1
src/tools/psb_cdalv.f90
Unescape Escape View File

@ -235,7 +235,6 @@ subroutine psb_cdalv(m, v, ictxt, desc_a, info, flag)
endif
! estimate local cols number
!!$ loc_col = int((psb_colrow_+1.d0)*loc_row)+1
loc_col = min(2*loc_row,m)
allocate(desc_a%loc_to_glob(loc_col),&
&desc_a%lprm(1),stat=info)

1
test/Fileread/df_sample.f90
Unescape Escape View File

@ -245,7 +245,6 @@ program df_sample
call psb_precset(pre,'asm',info,iv=(/novr,halo_,none_/))
case(rash_)
call psb_precset(pre,'asm',info,iv=(/novr,nohalo_,none_/))
case default
call psb_precset(pre,'ilu',info)
end select

Write Preview
Loading…
Cancel
Save