fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-06 22:55:12 +00:00
Code Issues Packages Projects Releases Wiki Activity

mld2p4:

Browse Source 
mlprec/mld_daggrmat_raw_asb.F90
mlprec/mld_dmlprec_aply.f90
mlprec/mld_dmlprec_bld.f90
mlprec/mld_prec_type.f90
mlprec/mld_zaggrmat_raw_asb.F90
mlprec/mld_zmlprec_aply.f90
mlprec/mld_zmlprec_bld.f90

Merged changes from mld2p4-intermesh branch up to rev. 2809.
This commit is contained in:
Salvatore Filippone
2008-02-08 12:36:53 +00:00
parent 5cd03fa531
commit cd1ed7f505
7 changed files with 178 additions and 364 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/mld_daggrmat_raw_asb.F90
+30
View File
@@ -99,6 +99,7 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
character(len=20) :: name
type(psb_dspmat_type) :: b
integer, pointer :: nzbr(:), idisp(:)
type(psb_dspmat_type), pointer :: am1,am2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzt, naggrm1, i
@@ -116,6 +117,13 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
am2 => p%av(mld_sm_pr_t_)
am1 => p%av(mld_sm_pr_)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
naggr = p%nlaggr(me+1)
ntaggr = sum(p%nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
@@ -140,6 +148,28 @@ subroutine mld_daggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
call psb_sp_all(ncol,ntaggr,am1,ncol,info)
else
call psb_sp_all(ncol,naggr,am1,ncol,info)
end if
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
end if
do i=1,nrow
am1%aspk(i) = done
am1%ia1(i) = i
am1%ia2(i) = p%mlia(i)
end do
am1%infoa(psb_nnz_) = nrow
call psb_spcnv(am1,info,afmt='csr',dupl=psb_dupl_add_)
call psb_transp(am1,am2)
call psb_sp_clip(a,b,info,jmax=nrow)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclip')
mlprec/mld_dmlprec_aply.f90
+51 -180
View File
@@ -436,8 +436,7 @@ contains
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%tx(max(n_row,n_col)),&
& mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info)
& stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/2*(nc2l+max(n_row,n_col)),0,0,0,0/),&
@@ -445,33 +444,15 @@ contains
goto 9999
end if
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(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(done,mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -510,22 +491,13 @@ contains
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_),mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(done,mlprec_wrk(ilev)%y2l,&
& done,mlprec_wrk(ilev-1)%y2l,&
& baseprecv(ilev)%map_desc,info,work=work)
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -704,9 +676,7 @@ contains
goto 9999
end if
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%x2l(:) = x
!
! STEP 2
!
@@ -758,29 +728,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%tx,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(done,mlprec_wrk(ilev-1)%tx,&
& dzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -826,24 +778,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%base_desc,info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,done,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(done,mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -1010,9 +951,6 @@ contains
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), 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)
@@ -1049,32 +987,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(done,mlprec_wrk(ilev-1)%x2l,&
& dzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -1135,25 +1052,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, dzero,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(done,mlprec_wrk(ilev+1)%y2l,&
& dzero,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -1358,10 +1263,6 @@ contains
& a_err='real(kind(1.d0))')
goto 9999
end if
mlprec_wrk(1)%x2l(:) = dzero
mlprec_wrk(1)%y2l(:) = dzero
mlprec_wrk(1)%tx(:) = dzero
mlprec_wrk(1)%ty(:) = dzero
call psb_geaxpby(done,x,dzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
@@ -1403,8 +1304,8 @@ contains
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%ty(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
if (info /= 0) then
info=4025
@@ -1413,30 +1314,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = dzero
mlprec_wrk(ilev)%y2l(:) = dzero
mlprec_wrk(ilev)%tx(:) = dzero
mlprec_wrk(ilev)%ty(:) = dzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%ty,dzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(done,mlprec_wrk(ilev-1)%ty,&
& dzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -1485,24 +1367,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(done,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, done,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(done,mlprec_wrk(ilev+1)%y2l,&
& done,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
mlprec/mld_dmlprec_bld.f90
+5
View File
@@ -153,6 +153,11 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info)
p%base_a => p%av(mld_ac_)
if (info==0) call psb_cdtransfer(desc_ac,p%desc_ac,info)
p%map_desc = psb_inter_desc(psb_map_aggr_,desc_a,&
& p%desc_ac,p%av(mld_sm_pr_t_),p%av(mld_sm_pr_))
! The two matrices from p%av() have been copied, may free them.
if (info == 0) call psb_sp_free(p%av(mld_sm_pr_t_),info)
if (info == 0) call psb_sp_free(p%av(mld_sm_pr_),info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdtransfer')
goto 9999
mlprec/mld_prec_type.f90
+5 -3
View File
@@ -61,7 +61,7 @@ module mld_prec_type
! blows up on some systems.
!
use psb_base_mod, only : psb_dspmat_type, psb_zspmat_type, psb_desc_type,&
& psb_sizeof
& psb_inter_desc_type, psb_sizeof
!
! Type: mld_dprec_type, mld_zprec_type
@@ -168,7 +168,7 @@ module mld_prec_type
type(psb_dspmat_type), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
real(kind(1.d0)), allocatable :: dorig(:)
type(psb_inter_desc_type) :: map_desc
end type mld_dbaseprc_type
type mld_dprec_type
@@ -187,7 +187,7 @@ module mld_prec_type
type(psb_zspmat_type), pointer :: base_a => null()
type(psb_desc_type), pointer :: base_desc => null()
complex(kind(1.d0)), allocatable :: dorig(:)
type(psb_inter_desc_type) :: map_desc
end type mld_zbaseprc_type
type mld_zprec_type
@@ -417,6 +417,7 @@ contains
val = val + psb_sizeof(prec%av(i))
end do
end if
val = val + psb_sizeof(prec%map_desc)
mld_dbaseprc_sizeof = val
@@ -453,6 +454,7 @@ contains
val = val + psb_sizeof(prec%av(i))
end do
end if
val = val + psb_sizeof(prec%map_desc)
mld_zbaseprc_sizeof = val
mlprec/mld_zaggrmat_raw_asb.F90
+30
View File
@@ -99,6 +99,7 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
character(len=20) :: name
type(psb_zspmat_type) :: b
integer, pointer :: nzbr(:), idisp(:)
type(psb_zspmat_type), pointer :: am1,am2
integer :: nrow, nglob, ncol, ntaggr, nzac, ip, ndx,&
& naggr, nzt,naggrm1, i
@@ -116,6 +117,13 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
nrow = psb_cd_get_local_rows(desc_a)
ncol = psb_cd_get_local_cols(desc_a)
am2 => p%av(mld_sm_pr_t_)
am1 => p%av(mld_sm_pr_)
call psb_nullify_sp(am1)
call psb_nullify_sp(am2)
naggr = p%nlaggr(me+1)
ntaggr = sum(p%nlaggr)
allocate(nzbr(np), idisp(np),stat=info)
@@ -140,6 +148,28 @@ subroutine mld_zaggrmat_raw_asb(a,desc_a,ac,desc_ac,p,info)
goto 9999
end if
if (p%iprcparm(mld_coarse_mat_) == mld_repl_mat_) then
call psb_sp_all(ncol,ntaggr,am1,ncol,info)
else
call psb_sp_all(ncol,naggr,am1,ncol,info)
end if
if (info /= 0) then
call psb_errpush(4010,name,a_err='spall')
goto 9999
end if
do i=1,nrow
am1%aspk(i) = zone
am1%ia1(i) = i
am1%ia2(i) = p%mlia(i)
end do
am1%infoa(psb_nnz_) = nrow
call psb_spcnv(am1,info,afmt='csr',dupl=psb_dupl_add_)
call psb_transp(am1,am2)
call psb_sp_clip(a,b,info,jmax=nrow)
if(info /= 0) then
call psb_errpush(4010,name,a_err='spclip')
mlprec/mld_zmlprec_aply.f90
+52 -181
View File
@@ -228,7 +228,7 @@ subroutine mld_zmlprec_aply(alpha,baseprecv,x,beta,y,desc_data,trans,work,info)
select case(baseprecv(2)%iprcparm(mld_smooth_pos_))
case(mld_post_smooth_)
select case (trans_)
case('N')
call mlt_post_ml_aply(alpha,baseprecv,x,beta,y,desc_data,trans_,work,info)
@@ -438,8 +438,7 @@ contains
nc2l = psb_cd_get_local_cols(baseprecv(ilev)%base_desc)
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
allocate(mlprec_wrk(ilev)%x2l(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%tx(max(n_row,n_col)),&
& mlprec_wrk(ilev)%ty(max(n_row,n_col)), stat=info)
& stat=info)
if (info /= 0) then
info=4025
call psb_errpush(info,name,i_err=(/2*(nc2l+max(n_row,n_col)),0,0,0,0/),&
@@ -447,33 +446,15 @@ contains
goto 9999
end if
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(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%x2l,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(zone,mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -512,22 +493,13 @@ contains
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_),mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(zone,mlprec_wrk(ilev)%y2l,&
& zone,mlprec_wrk(ilev-1)%y2l,&
& baseprecv(ilev)%map_desc,info,work=work)
end if
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -707,9 +679,7 @@ contains
goto 9999
end if
mlprec_wrk(1)%y2l(:) = zzero
mlprec_wrk(1)%x2l(:) = x
!
! STEP 2
!
@@ -761,29 +731,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%tx,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%tx,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(zone,mlprec_wrk(ilev-1)%tx,&
& zzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -829,24 +781,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,&
& baseprecv(ilev+1)%base_desc,info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(zone,mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -1014,9 +955,6 @@ contains
allocate(mlprec_wrk(1)%x2l(nc2l),mlprec_wrk(1)%y2l(nc2l), &
& mlprec_wrk(1)%tx(nc2l), 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)
@@ -1053,31 +991,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
if (debug_level >= psb_debug_inner_) &
& write(debug_unit,*) me,' ',trim(name), ' up sweep ', ilev
call psb_halo(mlprec_wrk(ilev-1)%x2l,&
& baseprecv(ilev-1)%base_desc,info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%x2l,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(zone,mlprec_wrk(ilev-1)%x2l,&
& zzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -1137,25 +1055,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l, zzero,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(zone,mlprec_wrk(ilev+1)%y2l,&
& zzero,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during prolongation')
goto 9999
@@ -1361,11 +1267,6 @@ contains
& a_err='real(kind(1.d0))')
goto 9999
end if
mlprec_wrk(1)%x2l(:) = zzero
mlprec_wrk(1)%y2l(:) = zzero
mlprec_wrk(1)%tx(:) = zzero
mlprec_wrk(1)%ty(:) = zzero
call psb_geaxpby(zone,x,zzero,mlprec_wrk(1)%x2l,&
& baseprecv(1)%base_desc,info)
@@ -1407,8 +1308,8 @@ contains
nr2l = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
ismth = baseprecv(ilev)%iprcparm(mld_aggr_kind_)
icm = baseprecv(ilev)%iprcparm(mld_coarse_mat_)
allocate(mlprec_wrk(ilev)%ty(nc2l),mlprec_wrk(ilev)%y2l(nc2l),&
& mlprec_wrk(ilev)%x2l(nc2l), stat=info)
allocate(mlprec_wrk(ilev)%tx(nc2l),mlprec_wrk(ilev)%ty(nc2l),&
& mlprec_wrk(ilev)%y2l(nc2l),mlprec_wrk(ilev)%x2l(nc2l), stat=info)
if (info /= 0) then
info=4025
@@ -1417,30 +1318,11 @@ contains
goto 9999
end if
mlprec_wrk(ilev)%x2l(:) = zzero
mlprec_wrk(ilev)%y2l(:) = zzero
mlprec_wrk(ilev)%tx(:) = zzero
mlprec_wrk(ilev)%ty(:) = zzero
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator transpose
!
call psb_halo(mlprec_wrk(ilev-1)%ty,baseprecv(ilev-1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev)%av(mld_sm_pr_t_),&
& mlprec_wrk(ilev-1)%ty,zzero,mlprec_wrk(ilev)%x2l,info)
else
!
! Apply the raw aggregation map transpose (take a shortcut)
!
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
! Apply prolongator transpose, i.e. restriction
call psb_forward_map(zone,mlprec_wrk(ilev-1)%ty,&
& zzero,mlprec_wrk(ilev)%x2l,&
& baseprecv(ilev)%map_desc,info,work=work)
if (info /=0) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
@@ -1489,24 +1371,13 @@ contains
ismth = baseprecv(ilev+1)%iprcparm(mld_aggr_kind_)
n_row = psb_cd_get_local_rows(baseprecv(ilev)%base_desc)
if (ismth /= mld_no_smooth_) then
!
! Apply the smoothed prolongator
!
if (ismth == mld_smooth_prol_) &
& call psb_halo(mlprec_wrk(ilev+1)%y2l,baseprecv(ilev+1)%base_desc,&
& info,work=work)
if (info == 0) call psb_csmm(zone,baseprecv(ilev+1)%av(mld_sm_pr_),&
& mlprec_wrk(ilev+1)%y2l,zone,mlprec_wrk(ilev)%y2l,info)
else
!
! Apply the raw aggregation map (take a shortcut)
!
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
!
! Apply prolongator
!
call psb_backward_map(zone,mlprec_wrk(ilev+1)%y2l,&
& zone,mlprec_wrk(ilev)%y2l,&
& baseprecv(ilev+1)%map_desc,info,work=work)
if (info /=0 ) then
call psb_errpush(4001,name,a_err='Error during restriction')
goto 9999
mlprec/mld_zmlprec_bld.f90
+5
View File
@@ -153,6 +153,11 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info)
p%base_a => p%av(mld_ac_)
if (info==0) call psb_cdtransfer(desc_ac,p%desc_ac,info)
p%map_desc = psb_inter_desc(psb_map_aggr_,desc_a,&
& p%desc_ac,p%av(mld_sm_pr_t_),p%av(mld_sm_pr_))
! The two matrices from p%av() have been copied, may free them.
if (info == 0) call psb_sp_free(p%av(mld_sm_pr_t_),info)
if (info == 0) call psb_sp_free(p%av(mld_sm_pr_),info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_cdtransfer')
goto 9999