Reworked preconditioner stuff. Currently contains some debug
statements. 



 base/modules/psb_c_mat_mod.f03
 base/modules/psb_d_base_mat_mod.f03
 base/modules/psb_d_mat_mod.f03
 base/modules/psb_error_mod.F90
 base/serial/f03/psb_d_csr_impl.f03
 base/serial/f03/psb_s_csr_impl.f03
 base/tools/psb_cspasb.f90
 base/tools/psb_dspasb.f90
 prec/Makefile
 prec/psb_c_bjacprec.f03
 prec/psb_c_diagprec.f03
 prec/psb_c_nullprec.f03
 prec/psb_cbjac_aply.f90
 prec/psb_cbjac_bld.f90
 prec/psb_cdiagsc_bld.f90
 prec/psb_cgprec_aply.f90
 prec/psb_cilu_fct.f90
 prec/psb_cprc_aply.f90
 prec/psb_cprecbld.f90
 prec/psb_cprecinit.f90
 prec/psb_cprecset.f90
 prec/psb_d_bjacprec.f03
 prec/psb_dbjac_aply.f90
 prec/psb_dbjac_bld.f90
 prec/psb_ddiagsc_bld.f90
 prec/psb_dgprec_aply.f90
 prec/psb_dprc_aply.f90
 prec/psb_dprecbld.f90
 prec/psb_dprecinit.f90
 prec/psb_dprecset.f90
 prec/psb_prec_type.f03
 prec/psb_s_bjacprec.f03
 prec/psb_s_diagprec.f03
 prec/psb_s_nullprec.f03
 prec/psb_sbjac_aply.f90
 prec/psb_sbjac_bld.f90
 prec/psb_sdiagsc_bld.f90
 prec/psb_sgprec_aply.f90
 prec/psb_silu_fct.f90
 prec/psb_sprc_aply.f90
 prec/psb_sprecbld.f90
 prec/psb_sprecinit.f90
 prec/psb_sprecset.f90
 prec/psb_z_bjacprec.f03
 prec/psb_z_diagprec.f03
 prec/psb_z_nullprec.f03
 prec/psb_zbjac_aply.f90
 prec/psb_zbjac_bld.f90
 prec/psb_zdiagsc_bld.f90
 prec/psb_zgprec_aply.f90
 prec/psb_zilu_fct.f90
 prec/psb_zprc_aply.f90
 prec/psb_zprecbld.f90
 prec/psb_zprecinit.f90
 prec/psb_zprecset.f90
 test/fileread/runs/cfs.inp
 test/fileread/runs/dfs.inp
 test/fileread/runs/zfs.inp
 test/pargen/Makefile
 test/pargen/runs/ppde.inp
 util/psb_mat_dist_mod.f90
psblas3-type-indexed
Salvatore Filippone 15 years ago
parent f3fd67a2ee
commit feeb610418

@ -1426,6 +1426,7 @@ contains
call a%set_dupl(psb_dupl_def_) call a%set_dupl(psb_dupl_def_)
end if end if
write(0,*)name,' ', present(mold), present(type),count( (/present(mold),present(type) /))
if (count( (/present(mold),present(type) /)) > 1) then if (count( (/present(mold),present(type) /)) > 1) then
info = 583 info = 583
call psb_errpush(info,name,a_err='TYPE, MOLD') call psb_errpush(info,name,a_err='TYPE, MOLD')

@ -1258,7 +1258,7 @@ contains
allocate(tmp(nac),stat=info) allocate(tmp(nac),stat=info)
if (info /= 0) info = 4000 if (info /= 0) info = 4000
if (info == 0) tmp(1:nac) = d(1:nac)*x(1:nac) if (info == 0) call inner_vscal(nac,d,x,tmp)
if (info == 0)& if (info == 0)&
& call a%base_cssm(alpha,tmp,beta,y,info,trans) & call a%base_cssm(alpha,tmp,beta,y,info,trans)
@ -1273,19 +1273,23 @@ contains
call psb_errpush(info,name,i_err=(/9,nar,0,0,0/)) call psb_errpush(info,name,i_err=(/9,nar,0,0,0/))
goto 9999 goto 9999
end if end if
allocate(tmp(nar),stat=info)
if (info /= 0) info = 4000
if (info == 0)&
& call a%base_cssm(done,x,dzero,tmp,info,trans)
if (info == 0) tmp(1:nar) = d(1:nar)*tmp(1:nar) if (beta == dzero) then
if (info == 0)& call a%base_cssm(alpha,x,dzero,y,info,trans)
& call psb_geaxpby(nar,alpha,tmp,beta,y,info) if (info == 0) call inner_vscal1(nar,d,y)
else
if (info == 0) then allocate(tmp(nar),stat=info)
deallocate(tmp,stat=info) if (info /= 0) info = 4000
if (info /= 0) info = 4000 if (info == 0)&
& call a%base_cssm(alpha,x,dzero,tmp,info,trans)
if (info == 0) call inner_vscal1(nar,d,tmp)
if (info == 0)&
& call psb_geaxpby(nar,done,tmp,beta,y,info)
if (info == 0) then
deallocate(tmp,stat=info)
if (info /= 0) info = 4000
end if
end if end if
else else
@ -1318,8 +1322,31 @@ contains
return return
end if end if
return return
contains
subroutine inner_vscal(n,d,x,y)
implicit none
integer, intent(in) :: n
real(psb_dpk_), intent(in) :: d(*),x(*)
real(psb_dpk_), intent(out) :: y(*)
integer :: i
do i=1,n
y(i) = d(i)*x(i)
end do
end subroutine inner_vscal
subroutine inner_vscal1(n,d,x)
implicit none
integer, intent(in) :: n
real(psb_dpk_), intent(in) :: d(*)
real(psb_dpk_), intent(inout) :: x(*)
integer :: i
do i=1,n
x(i) = d(i)*x(i)
end do
end subroutine inner_vscal1
end subroutine d_cssv end subroutine d_cssv

@ -1485,6 +1485,7 @@ contains
call a%set_dupl(psb_dupl_def_) call a%set_dupl(psb_dupl_def_)
end if end if
write(0,*)name,' ', present(mold), present(type),count( (/present(mold),present(type) /))
if (count( (/present(mold),present(type) /)) > 1) then if (count( (/present(mold),present(type) /)) > 1) then
info = 583 info = 583
call psb_errpush(info,name,a_err='TYPE, MOLD') call psb_errpush(info,name,a_err='TYPE, MOLD')

@ -490,7 +490,9 @@ contains
write (0,'("Invalid state for communication descriptor")') write (0,'("Invalid state for communication descriptor")')
case (1123) case (1123)
write (0,'("Invalid combined state for A and DESC_A")') write (0,'("Invalid combined state for A and DESC_A")')
case(1124:1999) case (1124)
write (0,'("Invalid state for object:",a)') trim(a_e_d)
case(1125:1999)
write (0,'("computational error. code: ",i0)')err_c write (0,'("computational error. code: ",i0)')err_c
case(2010) case(2010)
write (0,'("BLACS error. Number of processes=-1")') write (0,'("BLACS error. Number of processes=-1")')

@ -117,9 +117,9 @@ contains
do i=1,m do i=1,m
acc = dzero acc = dzero
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = -acc y(i) = acc
end do end do
else else
@ -149,11 +149,11 @@ contains
else if (alpha == -done) then else if (alpha == -done) then
do i=1,m do i=1,m
acc = dzero acc = y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = y(i) -acc y(i) = acc
end do end do
else else
@ -172,21 +172,21 @@ contains
if (alpha == done) then if (alpha == done) then
do i=1,m do i=1,m
acc = dzero acc = -y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc + val(j) * x(ja(j))
enddo enddo
y(i) = -y(i) + acc y(i) = acc
end do end do
else if (alpha == -done) then else if (alpha == -done) then
do i=1,m do i=1,m
acc = dzero acc = y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = -y(i) -acc y(i) = acc
end do end do
else else
@ -408,9 +408,9 @@ contains
do i=1,m do i=1,m
acc(1:nc) = dzero acc(1:nc) = dzero
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) - val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = -acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else else
@ -430,21 +430,21 @@ contains
if (alpha == done) then if (alpha == done) then
do i=1,m do i=1,m
acc(1:nc) = dzero acc(1:nc) = y(i,1:nc)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = y(i,1:nc) + acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else if (alpha == -done) then else if (alpha == -done) then
do i=1,m do i=1,m
acc(1:nc) = dzero acc(1:nc) = y(i,1:nc)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) - val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = y(i,1:nc) -acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else else
@ -629,6 +629,17 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
goto 9999 goto 9999
end if end if
if (size(x)<m) then
info = 36
call psb_errpush(info,name,i_err=(/3,m,0,0,0/))
goto 9999
end if
if (size(y)<m) then
info = 36
call psb_errpush(info,name,i_err=(/5,m,0,0,0/))
goto 9999
end if
if (alpha == dzero) then if (alpha == dzero) then
if (beta == dzero) then if (beta == dzero) then
@ -666,9 +677,9 @@ subroutine d_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,x,tmp) & a%irp,a%ja,a%val,x,tmp)
do i = 1, m
y(i) = alpha*tmp(i) + beta*y(i) call psb_geaxpby(m,alpha,tmp,beta,y,info)
end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)

@ -117,9 +117,9 @@ contains
do i=1,m do i=1,m
acc = szero acc = szero
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = -acc y(i) = acc
end do end do
else else
@ -149,11 +149,11 @@ contains
else if (alpha == -sone) then else if (alpha == -sone) then
do i=1,m do i=1,m
acc = szero acc = y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = y(i) -acc y(i) = acc
end do end do
else else
@ -172,21 +172,21 @@ contains
if (alpha == sone) then if (alpha == sone) then
do i=1,m do i=1,m
acc = szero acc = -y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc + val(j) * x(ja(j))
enddo enddo
y(i) = -y(i) + acc y(i) = acc
end do end do
else if (alpha == -sone) then else if (alpha == -sone) then
do i=1,m do i=1,m
acc = szero acc = y(i)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j)) acc = acc - val(j) * x(ja(j))
enddo enddo
y(i) = -y(i) -acc y(i) = acc
end do end do
else else
@ -408,9 +408,9 @@ contains
do i=1,m do i=1,m
acc(1:nc) = szero acc(1:nc) = szero
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) - val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = -acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else else
@ -430,21 +430,21 @@ contains
if (alpha == sone) then if (alpha == sone) then
do i=1,m do i=1,m
acc(1:nc) = szero acc(1:nc) = y(i,1:nc)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = y(i,1:nc) + acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else if (alpha == -sone) then else if (alpha == -sone) then
do i=1,m do i=1,m
acc(1:nc) = szero acc(1:nc) = y(i,1:nc)
do j=irp(i), irp(i+1)-1 do j=irp(i), irp(i+1)-1
acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) acc(1:nc) = acc(1:nc) - val(j) * x(ja(j),1:nc)
enddo enddo
y(i,1:nc) = y(i,1:nc) -acc(1:nc) y(i,1:nc) = acc(1:nc)
end do end do
else else
@ -621,10 +621,8 @@ subroutine s_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
endif endif
tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
m = a%get_nrows() m = a%get_nrows()
if (.not. (a%is_triangle())) then if (.not. (a%is_triangle())) then
info = 1121 info = 1121
call psb_errpush(info,name) call psb_errpush(info,name)
@ -676,11 +674,12 @@ subroutine s_csr_cssv_impl(alpha,a,x,beta,y,info,trans)
if (info /= 0) then if (info /= 0) then
return return
end if end if
call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
& a%irp,a%ja,a%val,x,tmp) & a%irp,a%ja,a%val,x,tmp)
do i = 1, m
y(i) = alpha*tmp(i) + beta*y(i) call psb_geaxpby(m,alpha,tmp,beta,y,info)
end do
end if end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
@ -840,6 +839,7 @@ subroutine s_csr_cssm_impl(alpha,a,x,beta,y,info,trans)
tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
m = a%get_nrows() m = a%get_nrows()
nc = min(size(x,2) , size(y,2)) nc = min(size(x,2) , size(y,2))
@ -1079,7 +1079,7 @@ subroutine s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
integer, intent(in), optional :: jmin,jmax, nzin integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale logical, intent(in), optional :: rscale,cscale
logical :: appens_, rscale_, cscale_ logical :: append_, rscale_, cscale_
integer :: nzin_, jmin_, jmax_, err_act, i integer :: nzin_, jmin_, jmax_, err_act, i
character(len=20) :: name='csget' character(len=20) :: name='csget'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -1104,11 +1104,11 @@ subroutine s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
end if end if
if (present(append)) then if (present(append)) then
appens_=append append_=append
else else
appens_=.false. append_=.false.
endif endif
if ((appens_).and.(present(nzin))) then if ((append_).and.(present(nzin))) then
nzin_ = nzin nzin_ = nzin
else else
nzin_ = 0 nzin_ = 0
@ -1129,7 +1129,7 @@ subroutine s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
goto 9999 goto 9999
end if end if
call csr_getptn(imin,imax,jmin_,jmax_,a,nz,ia,ja,nzin_,appens_,info,iren) call csr_getptn(imin,imax,jmin_,jmax_,a,nz,ia,ja,nzin_,append_,info,iren)
if (rscale_) then if (rscale_) then
do i=nzin_+1, nzin_+nz do i=nzin_+1, nzin_+nz
@ -1177,7 +1177,7 @@ contains
integer, optional :: iren(:) integer, optional :: iren(:)
integer :: nzin_, nza, idx,i,j,k, nzt, irw, lrw integer :: nzin_, nza, idx,i,j,k, nzt, irw, lrw
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
character(len=20) :: name='coo_getrow' character(len=20) :: name='csr_getptn'
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
@ -1255,7 +1255,7 @@ subroutine s_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
integer, intent(in), optional :: jmin,jmax, nzin integer, intent(in), optional :: jmin,jmax, nzin
logical, intent(in), optional :: rscale,cscale logical, intent(in), optional :: rscale,cscale
logical :: appens_, rscale_, cscale_ logical :: append_, rscale_, cscale_
integer :: nzin_, jmin_, jmax_, err_act, i integer :: nzin_, jmin_, jmax_, err_act, i
character(len=20) :: name='csget' character(len=20) :: name='csget'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
@ -1280,11 +1280,11 @@ subroutine s_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
end if end if
if (present(append)) then if (present(append)) then
appens_=append append_=append
else else
appens_=.false. append_=.false.
endif endif
if ((appens_).and.(present(nzin))) then if ((append_).and.(present(nzin))) then
nzin_ = nzin nzin_ = nzin
else else
nzin_ = 0 nzin_ = 0
@ -1305,7 +1305,7 @@ subroutine s_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
goto 9999 goto 9999
end if end if
call csr_getrow(imin,imax,jmin_,jmax_,a,nz,ia,ja,val,nzin_,appens_,info,& call csr_getrow(imin,imax,jmin_,jmax_,a,nz,ia,ja,val,nzin_,append_,info,&
& iren) & iren)
if (rscale_) then if (rscale_) then

@ -119,6 +119,7 @@ subroutine psb_cspasb(a,desc_a, info, afmt, upd, dupl,mold)
call a%set_ncols(n_col) call a%set_ncols(n_col)
end if end if
write(0,*)name,' ', present(mold), present(afmt),count( (/present(mold),present(afmt) /))
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)

@ -118,6 +118,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold)
call a%set_ncols(n_col) call a%set_ncols(n_col)
end if end if
write(0,*)name,' ', present(mold), present(afmt),count( (/present(mold),present(afmt) /))
call a%cscnv(info,type=afmt,dupl=dupl, mold=mold) call a%cscnv(info,type=afmt,dupl=dupl, mold=mold)

@ -2,22 +2,20 @@ include ../Make.inc
LIBDIR=../lib LIBDIR=../lib
HERE=. HERE=.
MODOBJS= psb_prec_type.o psb_prec_mod.o psb_d_diagprec.o psb_d_nullprec.o \ MODOBJS= psb_prec_type.o psb_prec_mod.o\
psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o psb_d_bjacprec.o \
psb_s_diagprec.o psb_s_nullprec.o psb_s_bjacprec.o \
psb_c_diagprec.o psb_c_nullprec.o psb_c_bjacprec.o \
psb_z_diagprec.o psb_z_nullprec.o psb_z_bjacprec.o
F90OBJS= psb_dilu_fct.o\ F90OBJS= psb_dilu_fct.o\
psb_dprecbld.o psb_dprecset.o psb_dprecinit.o \ psb_dprecbld.o psb_dprecset.o psb_dprecinit.o \
psb_sbjac_bld.o psb_silu_fct.o\ psb_silu_fct.o\
psb_sprecbld.o psb_sprecset.o psb_sprecinit.o \ psb_sprecbld.o psb_sprecset.o psb_sprecinit.o \
psb_sdiagsc_bld.o \ psb_cilu_fct.o\
psb_sprc_aply.o psb_sgprec_aply.o psb_sbjac_aply.o\
psb_cbjac_bld.o psb_cilu_fct.o\
psb_cprecbld.o psb_cprecset.o psb_cprecinit.o \ psb_cprecbld.o psb_cprecset.o psb_cprecinit.o \
psb_cdiagsc_bld.o \ psb_zilu_fct.o\
psb_cprc_aply.o psb_cgprec_aply.o psb_cbjac_aply.o\ psb_zprecbld.o psb_zprecset.o psb_zprecinit.o
psb_zbjac_bld.o psb_zilu_fct.o\
psb_zprecbld.o psb_zprecset.o psb_zprecinit.o \
psb_zdiagsc_bld.o \
psb_zprc_aply.o psb_zgprec_aply.o psb_zbjac_aply.o
LIBMOD=psb_prec_mod$(.mod) LIBMOD=psb_prec_mod$(.mod)
LOCAL_MODS=$(MODOBJS:.o=$(.mod)) LOCAL_MODS=$(MODOBJS:.o=$(.mod))
@ -37,7 +35,10 @@ $(OBJS): $(LIBDIR)/psb_base_mod$(.mod)
$(F90OBJS): $(MODOBJS) $(F90OBJS): $(MODOBJS)
psb_prec_mod.o: psb_prec_type.o psb_prec_mod.o: psb_prec_type.o
psb_s_bjacprec.o psb_s_diagprec.o psb_s_nullprec.o: psb_prec_type.o psb_prec_mod.o
psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_type.o psb_prec_mod.o psb_d_bjacprec.o psb_d_diagprec.o psb_d_nullprec.o: psb_prec_type.o psb_prec_mod.o
psb_c_bjacprec.o psb_c_diagprec.o psb_c_nullprec.o: psb_prec_type.o psb_prec_mod.o
psb_z_bjacprec.o psb_z_diagprec.o psb_z_nullprec.o: psb_prec_type.o psb_prec_mod.o
veryclean: clean veryclean: clean
/bin/rm -f $(LIBNAME) /bin/rm -f $(LIBNAME)

@ -0,0 +1,569 @@
module psb_c_bjacprec
use psb_prec_type
type, extends(psb_c_base_prec_type) :: psb_c_bjac_prec_type
integer, allocatable :: iprcparm(:)
type(psb_c_sparse_mat), allocatable :: av(:)
complex(psb_spk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => c_bjac_apply
procedure, pass(prec) :: precbld => c_bjac_precbld
procedure, pass(prec) :: precinit => c_bjac_precinit
procedure, pass(prec) :: c_base_precseti => c_bjac_precseti
procedure, pass(prec) :: c_base_precsetr => c_bjac_precsetr
procedure, pass(prec) :: c_base_precsetc => c_bjac_precsetc
procedure, pass(prec) :: precfree => c_bjac_precfree
procedure, pass(prec) :: precdescr => c_bjac_precdescr
procedure, pass(prec) :: sizeof => c_bjac_sizeof
end type psb_c_bjac_prec_type
character(len=15), parameter, private :: &
& fact_names(0:2)=(/'None ','ILU(n) ',&
& 'ILU(eps) '/)
contains
subroutine c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_c_bjac_prec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: alpha,beta
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
! Local variables
integer :: n_row,n_col
complex(psb_spk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name='c_bjac_prec_apply'
character(len=20) :: ch_err
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row = psb_cd_get_local_rows(desc_data)
n_col = psb_cd_get_local_cols(desc_data)
if (size(x) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/))
goto 9999
end if
if (size(y) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=conjg(prec%d),choice=psb_none_, work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
end select
if (info /=0) then
ch_err="psb_spsm"
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_apply
subroutine c_bjac_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_bjac_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_realloc(psb_ifpsz,prec%iprcparm,info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
end if
prec%iprcparm(:) = 0
prec%iprcparm(psb_p_type_) = psb_bjac_
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
prec%iprcparm(psb_ilu_fill_in_) = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precinit
subroutine c_bjac_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
use psb_prec_mod
Implicit None
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_c_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_c_csr_sparse_mat), allocatable :: lf, uf
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name='c_bjac_precbld'
character(len=20) :: ch_err
if(psb_get_errstatus() /= 0) return
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(prec%av)) then
if (size(prec%av) < psb_bp_ilu_avsz) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
endif
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = nrow_a
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(prec%d)) then
if (size(prec%d) < n_row) then
deallocate(prec%d)
endif
endif
if (.not.allocated(prec%d)) then
allocate(prec%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,prec%d,info)
if(info==0) then
call prec%av(psb_l_pr_)%mv_from(lf)
call prec%av(psb_u_pr_)%mv_from(uf)
call prec%av(psb_l_pr_)%set_asb()
call prec%av(psb_u_pr_)%set_asb()
call prec%av(psb_l_pr_)%trim()
call prec%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!!$ call prec%av(psb_l_pr_)%print(30+me)
!!$ call prec%av(psb_u_pr_)%print(40+me)
!!$ do i=1,n_row
!!$ write(50+me,*) i,prec%d(i)
!!$ end do
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precbld
subroutine c_bjac_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_bjac_precset'
call psb_erractionsave(err_act)
info = 0
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
select case(what)
case (psb_f_type_)
if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.(prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precseti
subroutine c_bjac_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precsetr
subroutine c_bjac_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precsetc
subroutine c_bjac_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, i
character(len=20) :: name='c_bjac_precfree'
call psb_erractionsave(err_act)
info = 0
if (allocated(prec%av)) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
end if
if (allocated(prec%d)) then
deallocate(prec%d,stat=info)
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precfree
subroutine c_bjac_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_c_bjac_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='c_bjac_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
write(iout_,*) 'Block Jacobi with: ',&
& fact_names(prec%iprcparm(psb_f_type_))
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_bjac_precdescr
function c_bjac_sizeof(prec) result(val)
use psb_base_mod
class(psb_c_bjac_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
if (allocated(prec%d)) then
val = val + 2*psb_sizeof_sp * size(prec%d)
endif
if (allocated(prec%av)) then
val = val + psb_sizeof(prec%av(psb_l_pr_))
val = val + psb_sizeof(prec%av(psb_u_pr_))
endif
return
end function c_bjac_sizeof
end module psb_c_bjacprec

@ -0,0 +1,375 @@
module psb_c_diagprec
use psb_prec_type
type, extends(psb_c_base_prec_type) :: psb_c_diag_prec_type
complex(psb_spk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => c_diag_apply
procedure, pass(prec) :: precbld => c_diag_precbld
procedure, pass(prec) :: precinit => c_diag_precinit
procedure, pass(prec) :: c_base_precseti => c_diag_precseti
procedure, pass(prec) :: c_base_precsetr => c_diag_precsetr
procedure, pass(prec) :: c_base_precsetc => c_diag_precsetc
procedure, pass(prec) :: precfree => c_diag_precfree
procedure, pass(prec) :: precdescr => c_diag_precdescr
procedure, pass(prec) :: sizeof => c_diag_sizeof
end type psb_c_diag_prec_type
contains
subroutine c_diag_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_c_diag_prec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(in) :: alpha, beta
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character :: trans_
character(len=20) :: name='c_diag_prec_apply'
complex(psb_spk_), pointer :: ww(:)
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the DIAG preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < nrow) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (present(trans)) then
trans_ = psb_toupper(trans)
else
trans_='N'
end if
select case(trans_)
case('N')
case('T','C')
case default
info=40
call psb_errpush(info,name,&
& i_err=(/6,0,0,0,0/),a_err=trans_)
goto 9999
end select
if (size(work) >= size(x)) then
ww => work
else
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4025,name,&
& i_err=(/size(x),0,0,0,0/),a_err='complex(psb_spk_)')
goto 9999
end if
end if
if (trans_=='C') then
ww(1:nrow) = x(1:nrow)*conjg(prec%d(1:nrow))
else
ww(1:nrow) = x(1:nrow)*prec%d(1:nrow)
endif
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
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_apply
subroutine c_diag_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_diag_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precinit
subroutine c_diag_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_c_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow,i
character(len=20) :: name='c_diag_precbld'
call psb_erractionsave(err_act)
info = 0
nrow = psb_cd_get_local_cols(desc_a)
if (allocated(prec%d)) then
if (size(prec%d) < nrow) then
deallocate(prec%d,stat=info)
end if
end if
if ((info == 0).and.(.not.allocated(prec%d))) then
allocate(prec%d(nrow), stat=info)
end if
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call a%get_diag(prec%d,info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='get_diag')
goto 9999
end if
do i=1,nrow
if (prec%d(i) == dzero) then
prec%d(i) = done
else
prec%d(i) = done/prec%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precbld
subroutine c_diag_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precseti
subroutine c_diag_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precsetr
subroutine c_diag_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precsetc
subroutine c_diag_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precfree
subroutine c_diag_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_c_diag_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='c_diag_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'Diagonal scaling'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_diag_precdescr
function c_diag_sizeof(prec) result(val)
use psb_base_mod
class(psb_c_diag_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
val = val + 2*psb_sizeof_sp * size(prec%d)
return
end function c_diag_sizeof
end module psb_c_diagprec

@ -0,0 +1,293 @@
module psb_c_nullprec
use psb_prec_type
type, extends(psb_c_base_prec_type) :: psb_c_null_prec_type
contains
procedure, pass(prec) :: apply => c_null_apply
procedure, pass(prec) :: precbld => c_null_precbld
procedure, pass(prec) :: precinit => c_null_precinit
procedure, pass(prec) :: c_base_precseti => c_null_precseti
procedure, pass(prec) :: c_base_precsetr => c_null_precsetr
procedure, pass(prec) :: c_base_precsetc => c_null_precsetc
procedure, pass(prec) :: precfree => c_null_precfree
procedure, pass(prec) :: precdescr => c_null_precdescr
procedure, pass(prec) :: sizeof => c_null_sizeof
end type psb_c_null_prec_type
contains
subroutine c_null_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_c_null_prec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(in) :: alpha, beta
complex(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='c_null_prec_apply'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
if (info /= 0 ) then
info = 4010
call psb_errpush(infoi,name,a_err="psb_geaxpby")
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_apply
subroutine c_null_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precinit
subroutine c_null_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_c_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow
character(len=20) :: name='c_null_precbld'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precbld
subroutine c_null_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precseti
subroutine c_null_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precsetr
subroutine c_null_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precsetc
subroutine c_null_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precfree
subroutine c_null_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_c_null_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='c_null_precset'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'No preconditioning'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine c_null_precdescr
function c_null_sizeof(prec) result(val)
use psb_base_mod
class(psb_c_null_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
return
end function c_null_sizeof
end module psb_c_nullprec

@ -1,166 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_cbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! 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.
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_cbjac_aply
implicit none
type(psb_desc_type), intent(in) :: desc_data
type(psb_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
complex(psb_spk_),intent(in) :: alpha,beta
character(len=1) :: trans
complex(psb_spk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,n_col
complex(psb_spk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row=desc_data%matrix_data(psb_n_row_)
n_col=desc_data%matrix_data(psb_n_col_)
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T')
call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
case('C')
call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=conjg(prec%d),choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cbjac_aply

@ -1,181 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_cbjac_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_cbjac_bld
implicit none
!
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. array Arguments ..
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_cprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_c_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
name='psb_cbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
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
select case(p%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(p%av)) then
if (size(p%av) < psb_bp_ilu_avsz) then
do i=1,size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = p%desc_data%matrix_data(psb_n_row_)
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,p%d,info)
if(info==0) then
call p%av(psb_l_pr_)%mv_from(lf)
call p%av(psb_u_pr_)%mv_from(uf)
call p%av(psb_l_pr_)%set_asb()
call p%av(psb_u_pr_)%set_asb()
call p%av(psb_l_pr_)%trim()
call p%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cbjac_bld

@ -1,120 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_cdiagsc_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_cdiagsc_bld
Implicit None
type(psb_c_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_cprec_type),intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
! Local scalars
Integer :: err, n_row, n_col,I,ictxt,&
& me,np,mglob,err_act
integer :: int_err(5)
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
call psb_info(ictxt, me, np)
! diagonal scaling
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_realloc')
goto 9999
end if
!
! Retrieve the diagonal entries of the matrix A
!
call a%get_diag(p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_getdiag'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!
! Copy into p%desc_data the descriptor associated to A
!
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
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
! it is set to one/a_ii
!
do i=1,n_row
if (p%d(i) == czero) then
p%d(i) = cone
else
p%d(i) = cone/p%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cdiagsc_bld

@ -1,140 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_cgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_cgprec_aply
implicit none
type(psb_desc_type),intent(in) :: desc_data
type(psb_cprec_type), intent(in) :: prec
complex(psb_spk_),intent(in) :: x(:)
complex(psb_spk_),intent(inout) :: y(:)
complex(psb_spk_),intent(in) :: alpha,beta
character(len=1) :: trans
complex(psb_spk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,int_err(5)
complex(psb_spk_), pointer :: ww(:)
character :: trans_
integer :: ictxt,np,me, err_act
character(len=20) :: name, ch_err
name='psb_cgprec_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
select case(prec%iprcparm(psb_p_type_))
case(psb_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(psb_diag_)
if (size(work) >= size(x)) then
ww => work
else
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_)
if (trans_=='C') then
ww(1:n_row) = x(1:n_row)*conjg(prec%d(1:n_row))
else
ww(1:n_row) = x(1:n_row)*prec%d(1:n_row)
endif
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(psb_bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans_,work,info)
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_cgprec_aply

@ -133,7 +133,7 @@ contains
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
call trw%allocate(0,0,info) call trw%allocate(0,0,1)
if(info /= 0) then if(info /= 0) then
info=4010 info=4010
ch_err='psb_sp_all' ch_err='psb_sp_all'

@ -50,7 +50,7 @@ subroutine psb_cprc_aply(prec,x,y,desc_data,info,trans, work)
integer :: ictxt,np,me,err_act integer :: ictxt,np,me,err_act
character(len=20) :: name character(len=20) :: name
name='psb_cprec_aply' name='psb_prc_aply'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -74,10 +74,12 @@ subroutine psb_cprc_aply(prec,x,y,desc_data,info,trans, work)
end if end if
call psb_gprec_aply(cone,prec,x,czero,y,desc_data,trans_,work_,info) if (.not.allocated(prec%prec)) then
info = 1124
! If the original distribution has an overlap we should fix that. call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call prec%prec%apply(cone,x,czero,y,desc_data,info,trans_,work=work_)
if (present(work)) then if (present(work)) then
else else
deallocate(work_) deallocate(work_)
@ -151,7 +153,7 @@ subroutine psb_cprc_aply1(prec,x,desc_data,info,trans)
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
complex(psb_spk_), pointer :: WW(:), w1(:) complex(psb_spk_), pointer :: WW(:), w1(:)
character(len=20) :: name character(len=20) :: name
name='psb_cprec1' name='psb_prc_aply1'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -164,12 +166,17 @@ subroutine psb_cprc_aply1(prec,x,desc_data,info,trans)
trans_='N' trans_='N'
end if end if
if (.not.allocated(prec%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
allocate(ww(size(x)),w1(size(x)),stat=info) allocate(ww(size(x)),w1(size(x)),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
call psb_cprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) call prec%prec%apply(cone,x,czero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999 if(info /=0) goto 9999
x(:) = ww(:) x(:) = ww(:)
deallocate(ww,W1) deallocate(ww,W1)

@ -82,51 +82,14 @@ subroutine psb_cprecbld(a,desc_a,p,info,upd)
! ALso should define symbolic names for the preconditioners. ! ALso should define symbolic names for the preconditioners.
! !
call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& if (.not.allocated(p%prec)) then
& psb_diag_,is_legal_prec) info = 1124
call psb_errpush(info,name,a_err="preconditioner")
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(psb_p_type_))
case (psb_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 (psb_diag_)
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (psb_bjac_)
call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
& psb_f_ilu_n_,is_legal_ml_fact)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
goto 9999
end if
case default
info=4010
ch_err='Unknown psb_p_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if
end select call p%prec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -33,6 +33,9 @@ subroutine psb_cprecinit(p,ptype,info)
use psb_base_mod use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_cprecinit use psb_prec_mod, psb_protect_name => psb_cprecinit
use psb_c_nullprec
use psb_c_diagprec
use psb_c_bjacprec
implicit none implicit none
type(psb_cprec_type), intent(inout) :: p type(psb_cprec_type), intent(inout) :: p
@ -40,35 +43,29 @@ subroutine psb_cprecinit(p,ptype,info)
integer, intent(out) :: info integer, intent(out) :: info
info = 0 info = 0
call psb_realloc(psb_ifpsz,p%iprcparm,info)
if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info)
if (info /= 0) return
p%iprcparm(:) = 0
if (allocated(p%prec) ) then
call p%prec%precfree(info)
if (info == 0) deallocate(p%prec,stat=info)
if (info /= 0) return
end if
select case(psb_toupper(ptype(1:len_trim(ptype)))) select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC') case ('NONE','NOPREC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_noprec_ allocate(psb_c_null_prec_type :: p%prec, stat=info)
p%iprcparm(psb_f_type_) = psb_f_none_
case ('DIAG') case ('DIAG')
p%iprcparm(:) = 0 allocate(psb_c_diag_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_diag_
p%iprcparm(psb_f_type_) = psb_f_none_
case ('BJAC') case ('BJAC')
p%iprcparm(:) = 0 allocate(psb_c_bjac_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_bjac_
p%iprcparm(psb_f_type_) = psb_f_ilu_n_
p%iprcparm(psb_ilu_fill_in_) = 0
case default case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"' write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2 info = 2
end select end select
if (info == 0) call p%prec%precinit(info)
end subroutine psb_cprecinit end subroutine psb_cprecinit

@ -37,30 +37,18 @@ subroutine psb_cprecseti(p,what,val,info)
type(psb_cprec_type), intent(inout) :: p type(psb_cprec_type), intent(inout) :: p
integer :: what, val integer :: what, val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0 info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
select case(what) call p%prec%precset(what,val,info)
case (psb_f_type_)
if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_cprecseti end subroutine psb_cprecseti
@ -75,32 +63,18 @@ subroutine psb_cprecsets(p,what,val,info)
integer :: what integer :: what
real(psb_spk_) :: val real(psb_spk_) :: val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
! call p%prec%precset(what,val,info)
! This will have to be changed if/when we put together an ILU(eps)
! factorization.
!
select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_cprecsets end subroutine psb_cprecsets

@ -222,7 +222,6 @@ contains
integer :: int_err(5) integer :: int_err(5)
character :: trans, unitd character :: trans, unitd
type(psb_d_csr_sparse_mat), allocatable :: lf, uf type(psb_d_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me integer :: ictxt,np,me
character(len=20) :: name='d_bjac_precbld' character(len=20) :: name='d_bjac_precbld'
@ -299,7 +298,6 @@ contains
end if end if
endif endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need ! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,prec%d,info) call psb_ilu_fct(a,lf,uf,prec%d,info)

@ -1,158 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! 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.
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dbjac_aply
implicit none
type(psb_desc_type), intent(in) :: desc_data
type(psb_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1) :: trans
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,n_col
real(psb_dpk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row=desc_data%matrix_data(psb_n_row_)
n_col=desc_data%matrix_data(psb_n_col_)
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dbjac_aply

@ -1,181 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dbjac_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dbjac_bld
implicit none
!
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. array Arguments ..
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_dprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_d_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
name='psb_dbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
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
select case(p%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(p%av)) then
if (size(p%av) < psb_bp_ilu_avsz) then
do i=1,size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = p%desc_data%matrix_data(psb_n_row_)
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,p%d,info)
if(info==0) then
call p%av(psb_l_pr_)%mv_from(lf)
call p%av(psb_u_pr_)%mv_from(uf)
call p%av(psb_l_pr_)%set_asb()
call p%av(psb_u_pr_)%set_asb()
call p%av(psb_l_pr_)%trim()
call p%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dbjac_bld

@ -1,120 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_ddiagsc_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_ddiagsc_bld
Implicit None
type(psb_d_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_dprec_type),intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
! Local scalars
Integer :: err, n_row, n_col,I,ictxt,&
& me,np,mglob, err_act
integer :: int_err(5)
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
call psb_info(ictxt, me, np)
! diagonal scaling
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_realloc')
goto 9999
end if
!
! Retrieve the diagonal entries of the matrix A
!
call a%get_diag(p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_getdiag'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!
! Copy into p%desc_data the descriptor associated to A
!
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
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
! it is set to one/a_ii
!
do i=1,n_row
if (p%d(i) == dzero) then
p%d(i) = done
else
p%d(i) = done/p%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_ddiagsc_bld

@ -1,135 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_dgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_dgprec_aply
implicit none
type(psb_desc_type),intent(in) :: desc_data
type(psb_dprec_type), intent(in) :: prec
real(psb_dpk_),intent(in) :: x(:)
real(psb_dpk_),intent(inout) :: y(:)
real(psb_dpk_),intent(in) :: alpha,beta
character(len=1) :: trans
real(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,int_err(5)
real(psb_dpk_), pointer :: ww(:)
character :: trans_
integer :: ictxt,np,me, err_act
character(len=20) :: name, ch_err
name='psb_dgprec_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
select case(prec%iprcparm(psb_p_type_))
case(psb_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(psb_diag_)
if (size(work) >= size(x)) then
ww => work
else
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(psb_bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans_,work,info)
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_dgprec_aply

@ -73,13 +73,12 @@ subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work)
end if end if
if (.not.allocated(prec%dprec)) then if (.not.allocated(prec%prec)) then
info = 1124 info = 1124
call psb_errpush(info,name,a_err="preconditioner") call psb_errpush(info,name,a_err="preconditioner")
goto 9999 goto 9999
end if end if
!!$ call psb_gprec_aply(done,prec,x,dzero,y,desc_data,trans_,work_,info) call prec%prec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_)
call prec%dprec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_)
if (present(work)) then if (present(work)) then
else else
deallocate(work_) deallocate(work_)
@ -153,7 +152,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
real(psb_dpk_), pointer :: WW(:), w1(:) real(psb_dpk_), pointer :: WW(:), w1(:)
character(len=20) :: name character(len=20) :: name
name='psb_dprec_aply1' name='psb_prc_aply1'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -166,7 +165,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
trans_='N' trans_='N'
end if end if
if (.not.allocated(prec%dprec)) then if (.not.allocated(prec%prec)) then
info = 1124 info = 1124
call psb_errpush(info,name,a_err="preconditioner") call psb_errpush(info,name,a_err="preconditioner")
goto 9999 goto 9999
@ -176,7 +175,7 @@ subroutine psb_dprc_aply1(prec,x,desc_data,info,trans)
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
call prec%dprec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1) call prec%prec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999 if(info /=0) goto 9999
x(:) = ww(:) x(:) = ww(:)
deallocate(ww,W1) deallocate(ww,W1)

@ -81,74 +81,26 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd)
! ALso should define symbolic names for the preconditioners. ! ALso should define symbolic names for the preconditioners.
! !
if (.false.) then if (.not.allocated(p%prec)) then
!!$ call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& info = 1124
!!$ & psb_diag_,is_legal_prec) call psb_errpush(info,name,a_err="preconditioner")
!!$ goto 9999
!!$ call psb_nullify_desc(p%desc_data)
!!$
!!$ select case(p%iprcparm(psb_p_type_))
!!$ case (psb_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 (psb_diag_)
!!$
!!$ call psb_diagsc_bld(a,desc_a,p,upd_,info)
!!$ if(info /= 0) then
!!$ info=4010
!!$ ch_err='psb_diagsc_bld'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$ end if
!!$
!!$ case (psb_bjac_)
!!$
!!$ call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
!!$ & psb_f_ilu_n_,is_legal_ml_fact)
!!$
!!$ call psb_bjac_bld(a,desc_a,p,upd_,info)
!!$
!!$ if(info /= 0) then
!!$ call psb_errpush(4010,name,a_err='psb_bjac_bld')
!!$ goto 9999
!!$ end if
!!$
!!$ case default
!!$ info=4010
!!$ ch_err='Unknown psb_p_type_'
!!$ call psb_errpush(info,name,a_err=ch_err)
!!$ goto 9999
!!$
!!$ end select
else
if (.not.allocated(p%dprec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call p%dprec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
end if end if
call p%prec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 continue 9999 continue
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then if (err_act == psb_act_abort_) then
call psb_error() call psb_error()
return
end if
return return
end if
return
end subroutine psb_dprecbld end subroutine psb_dprecbld

@ -43,57 +43,28 @@ subroutine psb_dprecinit(p,ptype,info)
info = 0 info = 0
if (.false.) then if (allocated(p%prec) ) then
!!$ call psb_realloc(psb_ifpsz,p%iprcparm,info) call p%prec%precfree(info)
!!$ if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info) if (info == 0) deallocate(p%prec,stat=info)
!!$ if (info /= 0) return if (info /= 0) return
!!$
!!$ select case(psb_toupper(ptype(1:len_trim(ptype))))
!!$ case ('NONE','NOPREC')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_noprec_
!!$ p%iprcparm(psb_f_type_) = psb_f_none_
!!$
!!$ case ('DIAG')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_diag_
!!$ p%iprcparm(psb_f_type_) = psb_f_none_
!!$
!!$ case ('BJAC')
!!$ p%iprcparm(:) = 0
!!$ p%iprcparm(psb_p_type_) = psb_bjac_
!!$ p%iprcparm(psb_f_type_) = psb_f_ilu_n_
!!$ p%iprcparm(psb_ilu_fill_in_) = 0
!!$
!!$ case default
!!$ write(0,*) 'Unknown preconditioner type request "',ptype,'"'
!!$ info = 2
!!$
!!$ end select
else
if (allocated(p%dprec) ) then
call p%dprec%precfree(info)
if (info == 0) deallocate(p%dprec,stat=info)
if (info /= 0) return
end if
select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
allocate(psb_d_null_prec_type :: p%dprec, stat=info)
case ('DIAG')
allocate(psb_d_diag_prec_type :: p%dprec, stat=info)
case ('BJAC')
allocate(psb_d_bjac_prec_type :: p%dprec, stat=info)
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
end select
if (info == 0) call p%dprec%precinit(info)
end if end if
select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC')
allocate(psb_d_null_prec_type :: p%prec, stat=info)
case ('DIAG')
allocate(psb_d_diag_prec_type :: p%prec, stat=info)
case ('BJAC')
allocate(psb_d_bjac_prec_type :: p%prec, stat=info)
case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2
end select
if (info == 0) call p%prec%precinit(info)
end subroutine psb_dprecinit end subroutine psb_dprecinit

@ -40,35 +40,15 @@ subroutine psb_dprecseti(p,what,val,info)
character(len=20) :: name='precset' character(len=20) :: name='precset'
info = 0 info = 0
if (.not.allocated(p%dprec)) then if (.not.allocated(p%prec)) then
info = 1124 info = 1124
call psb_errpush(info,name,a_err="preconditioner") call psb_errpush(info,name,a_err="preconditioner")
return return
!!$ goto 9999 !!$ goto 9999
end if end if
call p%dprec%precset(what,val,info) call p%prec%precset(what,val,info)
!!$ select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
!!$
!!$ case default
!!$ write(0,*) 'WHAT is invalid, ignoring user specification'
!!$
!!$ end select
return return
end subroutine psb_dprecseti end subroutine psb_dprecseti
@ -86,39 +66,15 @@ subroutine psb_dprecsetd(p,what,val,info)
character(len=20) :: name='precset' character(len=20) :: name='precset'
info = 0 info = 0
if (.not.allocated(p%dprec)) then if (.not.allocated(p%prec)) then
info = 1124 info = 1124
call psb_errpush(info,name,a_err="preconditioner") call psb_errpush(info,name,a_err="preconditioner")
return return
!!$ goto 9999 !!$ goto 9999
end if end if
call p%dprec%precset(what,val,info) call p%prec%precset(what,val,info)
!!$!
!!$! This will have to be changed if/when we put together an ILU(eps)
!!$! factorization.
!!$!
!!$ select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
!!$
!!$ case default
!!$ write(0,*) 'WHAT is invalid, ignoring user specification'
!!$
!!$ end select
return return
end subroutine psb_dprecsetd end subroutine psb_dprecsetd

File diff suppressed because it is too large Load Diff

@ -0,0 +1,563 @@
module psb_s_bjacprec
use psb_prec_type
type, extends(psb_s_base_prec_type) :: psb_s_bjac_prec_type
integer, allocatable :: iprcparm(:)
type(psb_s_sparse_mat), allocatable :: av(:)
real(psb_spk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => s_bjac_apply
procedure, pass(prec) :: precbld => s_bjac_precbld
procedure, pass(prec) :: precinit => s_bjac_precinit
procedure, pass(prec) :: s_base_precseti => s_bjac_precseti
procedure, pass(prec) :: s_base_precsetr => s_bjac_precsetr
procedure, pass(prec) :: s_base_precsetc => s_bjac_precsetc
procedure, pass(prec) :: precfree => s_bjac_precfree
procedure, pass(prec) :: precdescr => s_bjac_precdescr
procedure, pass(prec) :: sizeof => s_bjac_sizeof
end type psb_s_bjac_prec_type
character(len=15), parameter, private :: &
& fact_names(0:2)=(/'None ','ILU(n) ',&
& 'ILU(eps) '/)
contains
subroutine s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_s_bjac_prec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: alpha,beta
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
! Local variables
integer :: n_row,n_col
real(psb_spk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name='s_bjac_prec_apply'
character(len=20) :: ch_err
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row = psb_cd_get_local_rows(desc_data)
n_col = psb_cd_get_local_cols(desc_data)
if (size(x) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/))
goto 9999
end if
if (size(y) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T','C')
call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
end select
if (info /=0) then
ch_err="psb_spsm"
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_apply
subroutine s_bjac_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_realloc(psb_ifpsz,prec%iprcparm,info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
end if
prec%iprcparm(:) = 0
prec%iprcparm(psb_p_type_) = psb_bjac_
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
prec%iprcparm(psb_ilu_fill_in_) = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precinit
subroutine s_bjac_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
use psb_prec_mod
Implicit None
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_s_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_s_csr_sparse_mat), allocatable :: lf, uf
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name='s_bjac_precbld'
character(len=20) :: ch_err
if(psb_get_errstatus() /= 0) return
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(prec%av)) then
if (size(prec%av) < psb_bp_ilu_avsz) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
endif
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = nrow_a
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(prec%d)) then
if (size(prec%d) < n_row) then
deallocate(prec%d)
endif
endif
if (.not.allocated(prec%d)) then
allocate(prec%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,prec%d,info)
if(info==0) then
call prec%av(psb_l_pr_)%mv_from(lf)
call prec%av(psb_u_pr_)%mv_from(uf)
call prec%av(psb_l_pr_)%set_asb()
call prec%av(psb_u_pr_)%set_asb()
call prec%av(psb_l_pr_)%trim()
call prec%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!!$ call prec%av(psb_l_pr_)%print(30+me)
!!$ call prec%av(psb_u_pr_)%print(40+me)
!!$ do i=1,n_row
!!$ write(50+me,*) i,prec%d(i)
!!$ end do
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precbld
subroutine s_bjac_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_bjac_precset'
call psb_erractionsave(err_act)
info = 0
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
select case(what)
case (psb_f_type_)
if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.(prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precseti
subroutine s_bjac_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precsetr
subroutine s_bjac_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precsetc
subroutine s_bjac_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, i
character(len=20) :: name='s_bjac_precfree'
call psb_erractionsave(err_act)
info = 0
if (allocated(prec%av)) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
end if
if (allocated(prec%d)) then
deallocate(prec%d,stat=info)
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precfree
subroutine s_bjac_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_s_bjac_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='s_bjac_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
write(iout_,*) 'Block Jacobi with: ',&
& fact_names(prec%iprcparm(psb_f_type_))
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_bjac_precdescr
function s_bjac_sizeof(prec) result(val)
use psb_base_mod
class(psb_s_bjac_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
if (allocated(prec%d)) then
val = val + psb_sizeof_sp * size(prec%d)
endif
if (allocated(prec%av)) then
val = val + psb_sizeof(prec%av(psb_l_pr_))
val = val + psb_sizeof(prec%av(psb_u_pr_))
endif
return
end function s_bjac_sizeof
end module psb_s_bjacprec

@ -0,0 +1,352 @@
module psb_s_diagprec
use psb_prec_type
type, extends(psb_s_base_prec_type) :: psb_s_diag_prec_type
real(psb_spk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => s_diag_apply
procedure, pass(prec) :: precbld => s_diag_precbld
procedure, pass(prec) :: precinit => s_diag_precinit
procedure, pass(prec) :: s_base_precseti => s_diag_precseti
procedure, pass(prec) :: s_base_precsetr => s_diag_precsetr
procedure, pass(prec) :: s_base_precsetc => s_diag_precsetc
procedure, pass(prec) :: precfree => s_diag_precfree
procedure, pass(prec) :: precdescr => s_diag_precdescr
procedure, pass(prec) :: sizeof => s_diag_sizeof
end type psb_s_diag_prec_type
contains
subroutine s_diag_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_s_diag_prec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(in) :: alpha, beta
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='s_diag_prec_apply'
real(psb_spk_), pointer :: ww(:)
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the DIAG preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < nrow) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(work) >= size(x)) then
ww => work
else
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4025,name,i_err=(/size(x),0,0,0,0/),a_err='real(psb_spk_)')
goto 9999
end if
end if
ww(1:nrow) = x(1:nrow)*prec%d(1:nrow)
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
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_apply
subroutine s_diag_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_diag_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precinit
subroutine s_diag_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_s_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow,i
character(len=20) :: name='s_diag_precbld'
call psb_erractionsave(err_act)
info = 0
nrow = psb_cd_get_local_cols(desc_a)
if (allocated(prec%d)) then
if (size(prec%d) < nrow) then
deallocate(prec%d,stat=info)
end if
end if
if ((info == 0).and.(.not.allocated(prec%d))) then
allocate(prec%d(nrow), stat=info)
end if
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call a%get_diag(prec%d,info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='get_diag')
goto 9999
end if
do i=1,nrow
if (prec%d(i) == dzero) then
prec%d(i) = done
else
prec%d(i) = done/prec%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precbld
subroutine s_diag_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precseti
subroutine s_diag_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precsetr
subroutine s_diag_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precsetc
subroutine s_diag_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precfree
subroutine s_diag_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_s_diag_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='s_diag_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'Diagonal scaling'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_diag_precdescr
function s_diag_sizeof(prec) result(val)
use psb_base_mod
class(psb_s_diag_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
val = val + psb_sizeof_sp * size(prec%d)
return
end function s_diag_sizeof
end module psb_s_diagprec

@ -0,0 +1,293 @@
module psb_s_nullprec
use psb_prec_type
type, extends(psb_s_base_prec_type) :: psb_s_null_prec_type
contains
procedure, pass(prec) :: apply => s_null_apply
procedure, pass(prec) :: precbld => s_null_precbld
procedure, pass(prec) :: precinit => s_null_precinit
procedure, pass(prec) :: s_base_precseti => s_null_precseti
procedure, pass(prec) :: s_base_precsetr => s_null_precsetr
procedure, pass(prec) :: s_base_precsetc => s_null_precsetc
procedure, pass(prec) :: precfree => s_null_precfree
procedure, pass(prec) :: precdescr => s_null_precdescr
procedure, pass(prec) :: sizeof => s_null_sizeof
end type psb_s_null_prec_type
contains
subroutine s_null_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_s_null_prec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(in) :: alpha, beta
real(psb_spk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
real(psb_spk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='s_null_prec_apply'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
if (info /= 0 ) then
info = 4010
call psb_errpush(infoi,name,a_err="psb_geaxpby")
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_apply
subroutine s_null_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precinit
subroutine s_null_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_s_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow
character(len=20) :: name='s_null_precbld'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precbld
subroutine s_null_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precseti
subroutine s_null_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_spk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precsetr
subroutine s_null_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precsetc
subroutine s_null_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='s_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precfree
subroutine s_null_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_s_null_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='s_null_precset'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'No preconditioning'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine s_null_precdescr
function s_null_sizeof(prec) result(val)
use psb_base_mod
class(psb_s_null_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
return
end function s_null_sizeof
end module psb_s_nullprec

@ -1,158 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_sbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! 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.
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_sbjac_aply
implicit none
type(psb_desc_type), intent(in) :: desc_data
type(psb_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
real(psb_spk_),intent(in) :: alpha,beta
character(len=1) :: trans
real(psb_spk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,n_col
real(psb_spk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row=desc_data%matrix_data(psb_n_row_)
n_col=desc_data%matrix_data(psb_n_col_)
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T','C')
call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_sbjac_aply

@ -1,181 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_sbjac_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_sbjac_bld
implicit none
!
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. array Arguments ..
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_sprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_s_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
name='psb_sbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
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
select case(p%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(p%av)) then
if (size(p%av) < psb_bp_ilu_avsz) then
do i=1,size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = p%desc_data%matrix_data(psb_n_row_)
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,p%d,info)
if(info==0) then
call p%av(psb_l_pr_)%mv_from(lf)
call p%av(psb_u_pr_)%mv_from(uf)
call p%av(psb_l_pr_)%set_asb()
call p%av(psb_u_pr_)%set_asb()
call p%av(psb_l_pr_)%trim()
call p%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_sbjac_bld

@ -1,120 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_sdiagsc_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_sdiagsc_bld
Implicit None
type(psb_s_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_sprec_type),intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
! Local scalars
Integer :: err, n_row, n_col,I,ictxt,&
& me,np,mglob, err_act
integer :: int_err(5)
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
call psb_info(ictxt, me, np)
! diagonal scaling
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_realloc')
goto 9999
end if
!
! Retrieve the diagonal entries of the matrix A
!
call a%get_diag(p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_getdiag'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!
! Copy into p%desc_data the descriptor associated to A
!
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
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
! it is set to one/a_ii
!
do i=1,n_row
if (p%d(i) == szero) then
p%d(i) = sone
else
p%d(i) = sone/p%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_sdiagsc_bld

@ -1,135 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_sgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_sgprec_aply
implicit none
type(psb_desc_type),intent(in) :: desc_data
type(psb_sprec_type), intent(in) :: prec
real(psb_spk_),intent(in) :: x(:)
real(psb_spk_),intent(inout) :: y(:)
real(psb_spk_),intent(in) :: alpha,beta
character(len=1) :: trans
real(psb_spk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,int_err(5)
real(psb_spk_), pointer :: ww(:)
character :: trans_
integer :: ictxt,np,me, err_act
character(len=20) :: name, ch_err
name='psb_sgprec_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
select case(prec%iprcparm(psb_p_type_))
case(psb_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(psb_diag_)
if (size(work) >= size(x)) then
ww => work
else
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(psb_bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans_,work,info)
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_sgprec_aply

@ -136,7 +136,7 @@ contains
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
call trw%allocate(0,0,info) call trw%allocate(0,0,1)
if(info /= 0) then if(info /= 0) then
info=4010 info=4010
ch_err='psb_sp_all' ch_err='psb_sp_all'

@ -73,8 +73,12 @@ subroutine psb_sprc_aply(prec,x,y,desc_data,info,trans, work)
end if end if
call psb_gprec_aply(sone,prec,x,szero,y,desc_data,trans_,work_,info) if (.not.allocated(prec%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call prec%prec%apply(sone,x,szero,y,desc_data,info,trans_,work=work_)
if (present(work)) then if (present(work)) then
else else
deallocate(work_) deallocate(work_)
@ -148,7 +152,7 @@ subroutine psb_sprc_aply1(prec,x,desc_data,info,trans)
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
real(psb_spk_), pointer :: WW(:), w1(:) real(psb_spk_), pointer :: WW(:), w1(:)
character(len=20) :: name character(len=20) :: name
name='psb_sprec_aply1' name='psb_prc_aply1'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -161,12 +165,17 @@ subroutine psb_sprc_aply1(prec,x,desc_data,info,trans)
trans_='N' trans_='N'
end if end if
if (.not.allocated(prec%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
allocate(ww(size(x)),w1(size(x)),stat=info) allocate(ww(size(x)),w1(size(x)),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
call psb_sprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) call prec%prec%apply(sone,x,szero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999 if(info /=0) goto 9999
x(:) = ww(:) x(:) = ww(:)
deallocate(ww,W1) deallocate(ww,W1)

@ -81,51 +81,14 @@ subroutine psb_sprecbld(a,desc_a,p,info,upd)
! ALso should define symbolic names for the preconditioners. ! ALso should define symbolic names for the preconditioners.
! !
call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& if (.not.allocated(p%prec)) then
& psb_diag_,is_legal_prec) info = 1124
call psb_errpush(info,name,a_err="preconditioner")
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(psb_p_type_))
case (psb_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 (psb_diag_)
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (psb_bjac_)
call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
& psb_f_ilu_n_,is_legal_ml_fact)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
goto 9999
end if
case default
info=4010
ch_err='Unknown psb_p_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if
end select call p%prec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -33,6 +33,9 @@ subroutine psb_sprecinit(p,ptype,info)
use psb_base_mod use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_sprecinit use psb_prec_mod, psb_protect_name => psb_sprecinit
use psb_s_nullprec
use psb_s_diagprec
use psb_s_bjacprec
implicit none implicit none
type(psb_sprec_type), intent(inout) :: p type(psb_sprec_type), intent(inout) :: p
character(len=*), intent(in) :: ptype character(len=*), intent(in) :: ptype
@ -40,32 +43,28 @@ subroutine psb_sprecinit(p,ptype,info)
info = 0 info = 0
call psb_realloc(psb_ifpsz,p%iprcparm,info) if (allocated(p%prec) ) then
if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info) call p%prec%precfree(info)
if (info /= 0) return if (info == 0) deallocate(p%prec,stat=info)
p%iprcparm(:) = 0 if (info /= 0) return
end if
select case(psb_toupper(ptype(1:len_trim(ptype)))) select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC') case ('NONE','NOPREC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_noprec_ allocate(psb_s_null_prec_type :: p%prec, stat=info)
p%iprcparm(psb_f_type_) = psb_f_none_
case ('DIAG') case ('DIAG')
p%iprcparm(:) = 0 allocate(psb_s_diag_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_diag_
p%iprcparm(psb_f_type_) = psb_f_none_
case ('BJAC') case ('BJAC')
p%iprcparm(:) = 0 allocate(psb_s_bjac_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_bjac_
p%iprcparm(psb_f_type_) = psb_f_ilu_n_
p%iprcparm(psb_ilu_fill_in_) = 0
case default case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"' write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2 info = 2
end select end select
if (info == 0) call p%prec%precinit(info)
end subroutine psb_sprecinit end subroutine psb_sprecinit

@ -37,30 +37,18 @@ subroutine psb_sprecseti(p,what,val,info)
type(psb_sprec_type), intent(inout) :: p type(psb_sprec_type), intent(inout) :: p
integer :: what, val integer :: what, val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0 info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
select case(what) call p%prec%precset(what,val,info)
case (psb_f_type_)
if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_sprecseti end subroutine psb_sprecseti
@ -75,32 +63,18 @@ subroutine psb_sprecsets(p,what,val,info)
integer :: what integer :: what
real(psb_spk_) :: val real(psb_spk_) :: val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
! call p%prec%precset(what,val,info)
! This will have to be changed if/when we put together an ILU(eps)
! factorization.
!
select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_sprecsets end subroutine psb_sprecsets

@ -0,0 +1,569 @@
module psb_z_bjacprec
use psb_prec_type
type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type
integer, allocatable :: iprcparm(:)
type(psb_z_sparse_mat), allocatable :: av(:)
complex(psb_dpk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => z_bjac_apply
procedure, pass(prec) :: precbld => z_bjac_precbld
procedure, pass(prec) :: precinit => z_bjac_precinit
procedure, pass(prec) :: z_base_precseti => z_bjac_precseti
procedure, pass(prec) :: z_base_precsetr => z_bjac_precsetr
procedure, pass(prec) :: z_base_precsetc => z_bjac_precsetc
procedure, pass(prec) :: precfree => z_bjac_precfree
procedure, pass(prec) :: precdescr => z_bjac_precdescr
procedure, pass(prec) :: sizeof => z_bjac_sizeof
end type psb_z_bjac_prec_type
character(len=15), parameter, private :: &
& fact_names(0:2)=(/'None ','ILU(n) ',&
& 'ILU(eps) '/)
contains
subroutine z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_z_bjac_prec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: alpha,beta
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
! Local variables
integer :: n_row,n_col
complex(psb_dpk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name='z_bjac_prec_apply'
character(len=20) :: ch_err
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row = psb_cd_get_local_rows(desc_data)
n_col = psb_cd_get_local_cols(desc_data)
if (size(x) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/2,n_row,0,0,0/))
goto 9999
end if
if (size(y) < n_row) then
info = 36
call psb_errpush(info,name,i_err=(/3,n_row,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < n_row) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
case('T')
call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
case('C')
call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=conjg(prec%d),choice=psb_none_, work=aux)
if(info ==0) call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
end select
if (info /=0) then
ch_err="psb_spsm"
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_apply
subroutine z_bjac_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_bjac_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_realloc(psb_ifpsz,prec%iprcparm,info)
if (info /= 0) then
info = 4000
call psb_Errpush(info,name)
goto 9999
end if
prec%iprcparm(:) = 0
prec%iprcparm(psb_p_type_) = psb_bjac_
prec%iprcparm(psb_f_type_) = psb_f_ilu_n_
prec%iprcparm(psb_ilu_fill_in_) = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precinit
subroutine z_bjac_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
use psb_prec_mod
Implicit None
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_z_bjac_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_z_csr_sparse_mat), allocatable :: lf, uf
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name='z_bjac_precbld'
character(len=20) :: ch_err
if(psb_get_errstatus() /= 0) return
info = 0
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(prec%av)) then
if (size(prec%av) < psb_bp_ilu_avsz) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
endif
end if
if (.not.allocated(prec%av)) then
allocate(prec%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = nrow_a
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(prec%d)) then
if (size(prec%d) < n_row) then
deallocate(prec%d)
endif
endif
if (.not.allocated(prec%d)) then
allocate(prec%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,prec%d,info)
if(info==0) then
call prec%av(psb_l_pr_)%mv_from(lf)
call prec%av(psb_u_pr_)%mv_from(uf)
call prec%av(psb_l_pr_)%set_asb()
call prec%av(psb_u_pr_)%set_asb()
call prec%av(psb_l_pr_)%trim()
call prec%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!!$ call prec%av(psb_l_pr_)%print(30+me)
!!$ call prec%av(psb_u_pr_)%print(40+me)
!!$ do i=1,n_row
!!$ write(50+me,*) i,prec%d(i)
!!$ end do
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precbld
subroutine z_bjac_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_bjac_precset'
call psb_erractionsave(err_act)
info = 0
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
select case(what)
case (psb_f_type_)
if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.(prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',prec%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
prec%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precseti
subroutine z_bjac_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precsetr
subroutine z_bjac_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='c_bjac_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precsetc
subroutine z_bjac_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, i
character(len=20) :: name='z_bjac_precfree'
call psb_erractionsave(err_act)
info = 0
if (allocated(prec%av)) then
do i=1,size(prec%av)
call prec%av(i)%free()
enddo
deallocate(prec%av,stat=info)
end if
if (allocated(prec%d)) then
deallocate(prec%d,stat=info)
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precfree
subroutine z_bjac_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_z_bjac_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='z_bjac_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
if (.not.allocated(prec%iprcparm)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
write(iout_,*) 'Block Jacobi with: ',&
& fact_names(prec%iprcparm(psb_f_type_))
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_bjac_precdescr
function z_bjac_sizeof(prec) result(val)
use psb_base_mod
class(psb_z_bjac_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
if (allocated(prec%d)) then
val = val + 2*psb_sizeof_dp * size(prec%d)
endif
if (allocated(prec%av)) then
val = val + psb_sizeof(prec%av(psb_l_pr_))
val = val + psb_sizeof(prec%av(psb_u_pr_))
endif
return
end function z_bjac_sizeof
end module psb_z_bjacprec

@ -0,0 +1,375 @@
module psb_z_diagprec
use psb_prec_type
type, extends(psb_z_base_prec_type) :: psb_z_diag_prec_type
complex(psb_dpk_), allocatable :: d(:)
contains
procedure, pass(prec) :: apply => z_diag_apply
procedure, pass(prec) :: precbld => z_diag_precbld
procedure, pass(prec) :: precinit => z_diag_precinit
procedure, pass(prec) :: z_base_precseti => z_diag_precseti
procedure, pass(prec) :: z_base_precsetr => z_diag_precsetr
procedure, pass(prec) :: z_base_precsetc => z_diag_precsetc
procedure, pass(prec) :: precfree => z_diag_precfree
procedure, pass(prec) :: precdescr => z_diag_precdescr
procedure, pass(prec) :: sizeof => z_diag_sizeof
end type psb_z_diag_prec_type
contains
subroutine z_diag_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_z_diag_prec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(in) :: alpha, beta
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character :: trans_
character(len=20) :: name='z_diag_prec_apply'
complex(psb_dpk_), pointer :: ww(:)
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the DIAG preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
if (.not.allocated(prec%d)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (size(prec%d) < nrow) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner: D")
goto 9999
end if
if (present(trans)) then
trans_ = psb_toupper(trans)
else
trans_='N'
end if
select case(trans_)
case('N')
case('T','C')
case default
info=40
call psb_errpush(info,name,&
& i_err=(/6,0,0,0,0/),a_err=trans_)
goto 9999
end select
if (size(work) >= size(x)) then
ww => work
else
allocate(ww(size(x)),stat=info)
if (info /= 0) then
call psb_errpush(4025,name,&
& i_err=(/size(x),0,0,0,0/),a_err='complex(psb_dpk_)')
goto 9999
end if
end if
if (trans_=='C') then
ww(1:nrow) = x(1:nrow)*conjg(prec%d(1:nrow))
else
ww(1:nrow) = x(1:nrow)*prec%d(1:nrow)
endif
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
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_apply
subroutine z_diag_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_diag_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precinit
subroutine z_diag_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_z_diag_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow,i
character(len=20) :: name='z_diag_precbld'
call psb_erractionsave(err_act)
info = 0
nrow = psb_cd_get_local_cols(desc_a)
if (allocated(prec%d)) then
if (size(prec%d) < nrow) then
deallocate(prec%d,stat=info)
end if
end if
if ((info == 0).and.(.not.allocated(prec%d))) then
allocate(prec%d(nrow), stat=info)
end if
if (info /= 0) then
info = 4000
call psb_errpush(info,name)
goto 9999
end if
call a%get_diag(prec%d,info)
if (info /= 0) then
info = 4010
call psb_errpush(info,name, a_err='get_diag')
goto 9999
end if
do i=1,nrow
if (prec%d(i) == dzero) then
prec%d(i) = done
else
prec%d(i) = done/prec%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precbld
subroutine z_diag_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precseti
subroutine z_diag_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precsetr
subroutine z_diag_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precsetc
subroutine z_diag_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_diag_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precfree
subroutine z_diag_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_z_diag_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='z_diag_precdescr'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'Diagonal scaling'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_diag_precdescr
function z_diag_sizeof(prec) result(val)
use psb_base_mod
class(psb_z_diag_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
val = val + 2*psb_sizeof_dp * size(prec%d)
return
end function z_diag_sizeof
end module psb_z_diagprec

@ -0,0 +1,293 @@
module psb_z_nullprec
use psb_prec_type
type, extends(psb_z_base_prec_type) :: psb_z_null_prec_type
contains
procedure, pass(prec) :: apply => z_null_apply
procedure, pass(prec) :: precbld => z_null_precbld
procedure, pass(prec) :: precinit => z_null_precinit
procedure, pass(prec) :: z_base_precseti => z_null_precseti
procedure, pass(prec) :: z_base_precsetr => z_null_precsetr
procedure, pass(prec) :: z_base_precsetc => z_null_precsetc
procedure, pass(prec) :: precfree => z_null_precfree
procedure, pass(prec) :: precdescr => z_null_precdescr
procedure, pass(prec) :: sizeof => z_null_sizeof
end type psb_z_null_prec_type
contains
subroutine z_null_apply(alpha,prec,x,beta,y,desc_data,info,trans,work)
use psb_base_mod
type(psb_desc_type),intent(in) :: desc_data
class(psb_z_null_prec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(in) :: alpha, beta
complex(psb_dpk_),intent(inout) :: y(:)
integer, intent(out) :: info
character(len=1), optional :: trans
complex(psb_dpk_),intent(inout), optional, target :: work(:)
Integer :: err_act, nrow
character(len=20) :: name='z_null_prec_apply'
call psb_erractionsave(err_act)
!
! This is the base version and we should throw an error.
! Or should it be the NULL preonditioner???
!
info = 0
nrow = psb_cd_get_local_rows(desc_data)
if (size(x) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/2,nrow,0,0,0/))
goto 9999
end if
if (size(y) < nrow) then
info = 36
call psb_errpush(info,name,i_err=(/3,nrow,0,0,0/))
goto 9999
end if
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
if (info /= 0 ) then
info = 4010
call psb_errpush(infoi,name,a_err="psb_geaxpby")
goto 9999
end if
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_apply
subroutine z_null_precinit(prec,info)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_null_precinit'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precinit
subroutine z_null_precbld(a,desc_a,prec,info,upd)
use psb_base_mod
Implicit None
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in), target :: desc_a
class(psb_z_null_prec_type),intent(inout) :: prec
integer, intent(out) :: info
character, intent(in), optional :: upd
Integer :: err_act, nrow
character(len=20) :: name='z_null_precbld'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precbld
subroutine z_null_precseti(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
integer, intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precseti
subroutine z_null_precsetr(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
real(psb_dpk_), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precsetr
subroutine z_null_precsetc(prec,what,val,info)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type),intent(inout) :: prec
integer, intent(in) :: what
character(len=*), intent(in) :: val
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precsetc
subroutine z_null_precfree(prec,info)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type), intent(inout) :: prec
integer, intent(out) :: info
Integer :: err_act, nrow
character(len=20) :: name='z_null_precset'
call psb_erractionsave(err_act)
info = 0
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precfree
subroutine z_null_precdescr(prec,iout)
use psb_base_mod
Implicit None
class(psb_z_null_prec_type), intent(in) :: prec
integer, intent(in), optional :: iout
Integer :: err_act, nrow, info
character(len=20) :: name='z_null_precset'
integer :: iout_
call psb_erractionsave(err_act)
info = 0
if (present(iout)) then
iout_ = iout
else
iout_ = 6
end if
write(iout_,*) 'No preconditioning'
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine z_null_precdescr
function z_null_sizeof(prec) result(val)
use psb_base_mod
class(psb_z_null_prec_type), intent(in) :: prec
integer(psb_long_int_k_) :: val
val = 0
return
end function z_null_sizeof
end module psb_z_nullprec

@ -1,166 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zbjac_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! 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.
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zbjac_aply
implicit none
type(psb_desc_type), intent(in) :: desc_data
type(psb_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
complex(psb_dpk_),intent(in) :: alpha,beta
character(len=1) :: trans
complex(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,n_col
complex(psb_dpk_), pointer :: ww(:), aux(:)
integer :: ictxt,np,me, err_act, int_err(5)
integer :: debug_level, debug_unit
character :: trans_
character(len=20) :: name, ch_err
name='psb_bjac_aply'
info = 0
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = psb_cd_get_context(desc_data)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N','T','C')
! Ok
case default
call psb_errpush(40,name)
goto 9999
end select
n_row=desc_data%matrix_data(psb_n_row_)
n_col=desc_data%matrix_data(psb_n_col_)
if (n_col <= size(work)) then
ww => work(1:n_col)
if ((4*n_col+n_col) <= size(work)) then
aux => work(n_col+1:)
else
allocate(aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
else
allocate(ww(n_col),aux(4*n_col),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
select case(prec%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
select case(trans_)
case('N')
call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_,work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_u_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_, work=aux)
if(info /=0) goto 9999
case('T')
call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=prec%d,choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
case('C')
call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,ww,desc_data,info,&
& trans=trans_,scale='L',diag=conjg(prec%d),choice=psb_none_, work=aux)
if(info /=0) goto 9999
call psb_spsm(alpha,prec%av(psb_l_pr_),ww,beta,y,desc_data,info,&
& trans=trans_,scale='U',choice=psb_none_,work=aux)
if(info /=0) goto 9999
end select
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid factorization')
goto 9999
end select
call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (n_col <= size(work)) then
if ((4*n_col+n_col) <= size(work)) then
else
deallocate(aux)
endif
else
deallocate(ww,aux)
endif
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zbjac_aply

@ -1,181 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zbjac_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zbjac_bld
implicit none
!
! .. Scalar Arguments ..
integer, intent(out) :: info
! .. array Arguments ..
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_zprec_type), intent(inout) :: p
type(psb_desc_type), intent(in) :: desc_a
character, intent(in) :: upd
! .. Local Scalars ..
integer :: i, m
integer :: int_err(5)
character :: trans, unitd
type(psb_z_csr_sparse_mat), allocatable :: lf, uf
real(psb_dpk_) :: t1,t2,t3,t4,t5,t6, t7, t8
integer nztota, err_act, n_row, nrow_a,n_col, nhalo
integer :: ictxt,np,me
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
name='psb_zbjac_bld'
call psb_erractionsave(err_act)
ictxt=psb_cd_get_context(desc_a)
call psb_info(ictxt, me, np)
m = a%get_nrows()
if (m < 0) then
info = 10
int_err(1) = 1
int_err(2) = m
call psb_errpush(info,name,i_err=int_err)
goto 9999
endif
trans = 'N'
unitd = 'U'
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
select case(p%iprcparm(psb_f_type_))
case(psb_f_ilu_n_)
if (allocated(p%av)) then
if (size(p%av) < psb_bp_ilu_avsz) then
do i=1,size(p%av)
call p%av(i)%free()
enddo
deallocate(p%av,stat=info)
endif
end if
if (.not.allocated(p%av)) then
allocate(p%av(psb_max_avsz),stat=info)
if (info /= 0) then
call psb_errpush(4000,name)
goto 9999
end if
endif
nrow_a = psb_cd_get_local_rows(desc_a)
nztota = a%get_nzeros()
n_col = psb_cd_get_local_cols(desc_a)
nhalo = n_col-nrow_a
n_row = p%desc_data%matrix_data(psb_n_row_)
allocate(lf,uf,stat=info)
if (info == 0) call lf%allocate(n_row,n_row,nztota)
if (info == 0) call uf%allocate(n_row,n_row,nztota)
if(info/=0) then
info=4010
ch_err='psb_sp_all'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (allocated(p%d)) then
if (size(p%d) < n_row) then
deallocate(p%d)
endif
endif
if (.not.allocated(p%d)) then
allocate(p%d(n_row),stat=info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate')
goto 9999
end if
endif
t3 = psb_wtime()
! This is where we have no renumbering, thus no need
call psb_ilu_fct(a,lf,uf,p%d,info)
if(info==0) then
call p%av(psb_l_pr_)%mv_from(lf)
call p%av(psb_u_pr_)%mv_from(uf)
call p%av(psb_l_pr_)%set_asb()
call p%av(psb_u_pr_)%set_asb()
call p%av(psb_l_pr_)%trim()
call p%av(psb_u_pr_)%trim()
else
info=4010
ch_err='psb_ilu_fct'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case(psb_f_none_)
info=4010
ch_err='Inconsistent prec psb_f_none_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
case default
info=4010
ch_err='Unknown psb_f_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zbjac_bld

@ -1,120 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zdiagsc_bld(a,desc_a,p,upd,info)
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zdiagsc_bld
Implicit None
type(psb_z_sparse_mat), intent(in), target :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_zprec_type),intent(inout) :: p
character, intent(in) :: upd
integer, intent(out) :: info
! Local scalars
Integer :: err, n_row, n_col,I,ictxt,&
& me,np,mglob,err_act
integer :: int_err(5)
integer,parameter :: iroot=psb_root_,iout=60,ilout=40
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=0
err=0
call psb_erractionsave(err_act)
name = 'psb_diagsc_bld'
info = 0
int_err(1) = 0
ictxt = psb_cd_get_context(desc_a)
n_row = psb_cd_get_local_rows(desc_a)
n_col = psb_cd_get_local_cols(desc_a)
mglob = psb_cd_get_global_rows(desc_a)
call psb_info(ictxt, me, np)
! diagonal scaling
call psb_realloc(n_col,p%d,info)
if (info /= 0) then
call psb_errpush(4010,name,a_err='psb_realloc')
goto 9999
end if
!
! Retrieve the diagonal entries of the matrix A
!
call a%get_diag(p%d,info)
if(info /= 0) then
info=4010
ch_err='psb_sp_getdiag'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
!
! Copy into p%desc_data the descriptor associated to A
!
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
!
! The i-th diagonal entry of the preconditioner is set to one if the
! corresponding entry a_ii of the sparse matrix A is zero; otherwise
! it is set to one/a_ii
!
do i=1,n_row
if (p%d(i) == zzero) then
p%d(i) = zone
else
p%d(i) = zone/p%d(i)
endif
end do
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zdiagsc_bld

@ -1,140 +0,0 @@
!!$
!!$ Parallel Sparse BLAS version 2.2
!!$ (C) Copyright 2006/2007/2008
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari University of Rome Tor Vergata
!!$
!!$ Redistribution and use in source and binary forms, with or without
!!$ modification, are permitted provided that the following conditions
!!$ are met:
!!$ 1. Redistributions of source code must retain the above copyright
!!$ notice, this list of conditions and the following disclaimer.
!!$ 2. Redistributions in binary form must reproduce the above copyright
!!$ notice, this list of conditions, and the following disclaimer in the
!!$ documentation and/or other materials provided with the distribution.
!!$ 3. The name of the PSBLAS group or the names of its contributors may
!!$ not be used to endorse or promote products derived from this
!!$ software without specific written permission.
!!$
!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!$ POSSIBILITY OF SUCH DAMAGE.
!!$
!!$
subroutine psb_zgprec_aply(alpha,prec,x,beta,y,desc_data,trans,work,info)
!
! Compute Y <- beta*Y + alpha*K^-1 X
! where K is a a basic preconditioner stored in prec
!
use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zgprec_aply
implicit none
type(psb_desc_type),intent(in) :: desc_data
type(psb_zprec_type), intent(in) :: prec
complex(psb_dpk_),intent(in) :: x(:)
complex(psb_dpk_),intent(inout) :: y(:)
complex(psb_dpk_),intent(in) :: alpha,beta
character(len=1) :: trans
complex(psb_dpk_),target :: work(:)
integer, intent(out) :: info
! Local variables
integer :: n_row,int_err(5)
complex(psb_dpk_), pointer :: ww(:)
character :: trans_
integer :: ictxt,np,me, err_act
character(len=20) :: name, ch_err
name='psb_zgprec_aply'
info = 0
call psb_erractionsave(err_act)
ictxt=desc_data%matrix_data(psb_ctxt_)
call psb_info(ictxt, me, np)
trans_ = psb_toupper(trans)
select case(trans_)
case('N')
case('T','C')
case default
info=40
int_err(1)=6
ch_err(2:2)=trans
goto 9999
end select
select case(prec%iprcparm(psb_p_type_))
case(psb_noprec_)
call psb_geaxpby(alpha,x,beta,y,desc_data,info)
case(psb_diag_)
if (size(work) >= size(x)) then
ww => work
else
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_)
if (trans_=='C') then
ww(1:n_row) = x(1:n_row)*conjg(prec%d(1:n_row))
else
ww(1:n_row) = x(1:n_row)*prec%d(1:n_row)
endif
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(psb_bjac_)
call psb_bjac_aply(alpha,prec,x,beta,y,desc_data,trans_,work,info)
if(info /= 0) then
info=4010
ch_err='psb_bjac_aply'
goto 9999
end if
case default
info = 4001
call psb_errpush(info,name,a_err='Invalid prectype')
goto 9999
end select
call psb_erractionrestore(err_act)
return
9999 continue
call psb_errpush(info,name,i_err=int_err,a_err=ch_err)
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_zgprec_aply

@ -133,7 +133,7 @@ contains
if(psb_get_errstatus() /= 0) return if(psb_get_errstatus() /= 0) return
info=0 info=0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
call trw%allocate(0,0,info) call trw%allocate(0,0,1)
if(info /= 0) then if(info /= 0) then
info=4010 info=4010
ch_err='psb_sp_all' ch_err='psb_sp_all'

@ -50,7 +50,7 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
integer :: ictxt,np,me,err_act integer :: ictxt,np,me,err_act
character(len=20) :: name character(len=20) :: name
name='psb_zprec_aply' name='psb_prc_aply'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -74,10 +74,12 @@ subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work)
end if end if
call psb_gprec_aply(zone,prec,x,zzero,y,desc_data,trans_,work_,info) if (.not.allocated(prec%prec)) then
info = 1124
! If the original distribution has an overlap we should fix that. call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
call prec%prec%apply(zone,x,zzero,y,desc_data,info,trans_,work=work_)
if (present(work)) then if (present(work)) then
else else
deallocate(work_) deallocate(work_)
@ -151,7 +153,7 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
integer :: ictxt,np,me, err_act integer :: ictxt,np,me, err_act
complex(psb_dpk_), pointer :: WW(:), w1(:) complex(psb_dpk_), pointer :: WW(:), w1(:)
character(len=20) :: name character(len=20) :: name
name='psb_zprec1' name='psb_prc_aply1'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
@ -164,12 +166,17 @@ subroutine psb_zprc_aply1(prec,x,desc_data,info,trans)
trans_='N' trans_='N'
end if end if
if (.not.allocated(prec%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
goto 9999
end if
allocate(ww(size(x)),w1(size(x)),stat=info) allocate(ww(size(x)),w1(size(x)),stat=info)
if (info /= 0) then if (info /= 0) then
call psb_errpush(4010,name,a_err='Allocate') call psb_errpush(4010,name,a_err='Allocate')
goto 9999 goto 9999
end if end if
call psb_zprc_aply(prec,x,ww,desc_data,info,trans_,work=w1) call prec%prec%apply(zone,x,zzero,ww,desc_data,info,trans_,work=w1)
if(info /=0) goto 9999 if(info /=0) goto 9999
x(:) = ww(:) x(:) = ww(:)
deallocate(ww,W1) deallocate(ww,W1)

@ -82,51 +82,14 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd)
! ALso should define symbolic names for the preconditioners. ! ALso should define symbolic names for the preconditioners.
! !
call psb_check_def(p%iprcparm(psb_p_type_),'base_prec',& if (.not.allocated(p%prec)) then
& psb_diag_,is_legal_prec) info = 1124
call psb_errpush(info,name,a_err="preconditioner")
call psb_nullify_desc(p%desc_data)
select case(p%iprcparm(psb_p_type_))
case (psb_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 (psb_diag_)
call psb_diagsc_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
info=4010
ch_err='psb_diagsc_bld'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
case (psb_bjac_)
call psb_check_def(p%iprcparm(psb_f_type_),'fact',&
& psb_f_ilu_n_,is_legal_ml_fact)
call psb_bjac_bld(a,desc_a,p,upd_,info)
if(info /= 0) then
call psb_errpush(4010,name,a_err='psb_bjac_bld')
goto 9999
end if
case default
info=4010
ch_err='Unknown psb_p_type_'
call psb_errpush(info,name,a_err=ch_err)
goto 9999 goto 9999
end if
end select call p%prec%precbld(a,desc_a,info,upd)
if (info /= 0) goto 9999
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -33,6 +33,9 @@ subroutine psb_zprecinit(p,ptype,info)
use psb_base_mod use psb_base_mod
use psb_prec_mod, psb_protect_name => psb_zprecinit use psb_prec_mod, psb_protect_name => psb_zprecinit
use psb_z_nullprec
use psb_z_diagprec
use psb_z_bjacprec
implicit none implicit none
type(psb_zprec_type), intent(inout) :: p type(psb_zprec_type), intent(inout) :: p
@ -40,35 +43,29 @@ subroutine psb_zprecinit(p,ptype,info)
integer, intent(out) :: info integer, intent(out) :: info
info = 0 info = 0
call psb_realloc(psb_ifpsz,p%iprcparm,info)
if (info == 0) call psb_realloc(psb_rfpsz,p%rprcparm,info)
if (info /= 0) return
p%iprcparm(:) = 0
if (allocated(p%prec) ) then
call p%prec%precfree(info)
if (info == 0) deallocate(p%prec,stat=info)
if (info /= 0) return
end if
select case(psb_toupper(ptype(1:len_trim(ptype)))) select case(psb_toupper(ptype(1:len_trim(ptype))))
case ('NONE','NOPREC') case ('NONE','NOPREC')
p%iprcparm(:) = 0
p%iprcparm(psb_p_type_) = psb_noprec_ allocate(psb_z_null_prec_type :: p%prec, stat=info)
p%iprcparm(psb_f_type_) = psb_f_none_
case ('DIAG') case ('DIAG')
p%iprcparm(:) = 0 allocate(psb_z_diag_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_diag_
p%iprcparm(psb_f_type_) = psb_f_none_
case ('BJAC') case ('BJAC')
p%iprcparm(:) = 0 allocate(psb_z_bjac_prec_type :: p%prec, stat=info)
p%iprcparm(psb_p_type_) = psb_bjac_
p%iprcparm(psb_f_type_) = psb_f_ilu_n_
p%iprcparm(psb_ilu_fill_in_) = 0
case default case default
write(0,*) 'Unknown preconditioner type request "',ptype,'"' write(0,*) 'Unknown preconditioner type request "',ptype,'"'
info = 2 info = 2
end select end select
if (info == 0) call p%prec%precinit(info)
end subroutine psb_zprecinit end subroutine psb_zprecinit

@ -37,30 +37,18 @@ subroutine psb_zprecseti(p,what,val,info)
type(psb_zprec_type), intent(inout) :: p type(psb_zprec_type), intent(inout) :: p
integer :: what, val integer :: what, val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0 info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
select case(what) call p%prec%precset(what,val,info)
case (psb_f_type_)
if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_f_type_) = val
case (psb_ilu_fill_in_)
if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
& 'ignoring user specification'
return
endif
p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_zprecseti end subroutine psb_zprecseti
@ -75,32 +63,18 @@ subroutine psb_zprecsetd(p,what,val,info)
integer :: what integer :: what
real(psb_dpk_) :: val real(psb_dpk_) :: val
integer, intent(out) :: info integer, intent(out) :: info
character(len=20) :: name='precset'
info = 0
if (.not.allocated(p%prec)) then
info = 1124
call psb_errpush(info,name,a_err="preconditioner")
return
!!$ goto 9999
end if
! call p%prec%precset(what,val,info)
! This will have to be changed if/when we put together an ILU(eps)
! factorization.
!
select case(what)
!!$ case (psb_f_type_)
!!$ if (p%iprcparm(psb_p_type_) /= psb_bjac_) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_f_type_) = val
!!$
!!$ case (psb_ilu_fill_in_)
!!$ if ((p%iprcparm(psb_p_type_) /= psb_bjac_).or.(p%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then
!!$ write(0,*) 'WHAT is invalid for current preconditioner ',p%iprcparm(psb_p_type_),&
!!$ & 'ignoring user specification'
!!$ return
!!$ endif
!!$ p%iprcparm(psb_ilu_fill_in_) = val
case default
write(0,*) 'WHAT is invalid, ignoring user specification'
end select
return return
end subroutine psb_zprecsetd end subroutine psb_zprecsetd

@ -1,9 +1,9 @@
11 Number of inputs 11 Number of inputs
young1c.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or waveguide3D.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or
NONE http://www.cise.ufl.edu/research/sparse/matrices/index.html NONE http://www.cise.ufl.edu/research/sparse/matrices/index.html
MM File format: MM: Matrix Market HB: Harwell-Boeing. MM File format: MM: Matrix Market HB: Harwell-Boeing.
BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG
BJAC Preconditioner NONE DIAG BJAC DIAG Preconditioner NONE DIAG BJAC
CSR Storage format CSR COO JAD CSR Storage format CSR COO JAD
0 IPART: Partition method 0: BLK 2: graph (with Metis) 0 IPART: Partition method 0: BLK 2: graph (with Metis)
2 ISTOPC 2 ISTOPC

@ -1,13 +1,13 @@
11 Number of inputs 11 Number of inputs
thm200x120.mtx sherman3.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or 10974x10974.mtx sherman3.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or
NONE sherman3_rhs1.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html NONE sherman3_rhs1.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html
MM File format: MM: Matrix Market HB: Harwell-Boeing. MM File format: MM: Matrix Market HB: Harwell-Boeing.
BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG
BJAC Preconditioner NONE DIAG BJAC BJAC Preconditioner NONE DIAG BJAC
CSR Storage format CSR COO JAD CSR Storage format CSR COO JAD
2 IPART: Partition method 0: BLK 2: graph (with Metis) 0 IPART: Partition method 0: BLK 2: graph (with Metis)
2 ISTOPC 2 ISTOPC
01000 ITMAX 00010 ITMAX
01 ITRACE 01 ITRACE
30 IRST (restart for RGMRES and BiCGSTABL) 30 IRST (restart for RGMRES and BiCGSTABL)
1.d-6 EPS 1.d-6 EPS

@ -5,7 +5,7 @@ MM File format: MM: Matrix Market HB: Harwell-Boeing.
BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG BICGSTAB Iterative method: BiCGSTAB CGS RGMRES BiCGSTABL BICG CG
BJAC Preconditioner NONE DIAG BJAC BJAC Preconditioner NONE DIAG BJAC
CSR Storage format CSR COO JAD CSR Storage format CSR COO JAD
0 IPART: Partition method 0: BLK 2: graph (with Metis) 2 IPART: Partition method 0: BLK 2: graph (with Metis)
2 ISTOPC 2 ISTOPC
05000 ITMAX 05000 ITMAX
01 ITRACE 01 ITRACE

@ -17,7 +17,7 @@ EXEDIR=./runs
all: ppde spde all: ppde spde
ppde: ppde.o psb_d_csc_impl.o psb_d_csc_mat_mod.o ppde: ppde.o psb_d_csc_impl.o psb_d_csc_mat_mod.o
$(F90LINK) ppde.o psb_d_csc_impl.o psb_d_csc_mat_mod.o -o ppde $(PSBLAS_LIB) $(LDLIBS) $(F90LINK) -pg ppde.o psb_d_csc_impl.o psb_d_csc_mat_mod.o -o ppde $(PSBLAS_LIB) $(LDLIBS)
/bin/mv ppde $(EXEDIR) /bin/mv ppde $(EXEDIR)
psb_d_csc_impl.o ppde.o: psb_d_csc_mat_mod.o psb_d_csc_impl.o ppde.o: psb_d_csc_mat_mod.o

@ -4,8 +4,8 @@ BJAC Preconditioner NONE DIAG BJAC
CSR Storage format for matrix A: CSR COO JAD CSR Storage format for matrix A: CSR COO JAD
060 Domain size (acutal system is this**3) 060 Domain size (acutal system is this**3)
2 Stopping criterion 2 Stopping criterion
0400 MAXIT 0100 MAXIT
001 ITRACE 01 ITRACE
20 IRST restart for RGMRES and BiCGSTABL 20 IRST restart for RGMRES and BiCGSTABL

@ -1324,6 +1324,7 @@ contains
goto 9999 goto 9999
end if end if
write(0,*)name,' Calling spasb',psb_dupl_err_,' ',afmt
call psb_barrier(ictxt) call psb_barrier(ictxt)
t2 = psb_wtime() t2 = psb_wtime()
call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt)

Loading…
Cancel
Save