psblas2-dev:


			
			
				psblas3-type-indexed
			
			
		
Salvatore Filippone 17 years ago
parent 486a011f4b
commit 93b9ebbcf0

@ -586,9 +586,116 @@ Module psb_tools_mod
module procedure psb_dlinmap_asb, psb_zlinmap_asb module procedure psb_dlinmap_asb, psb_zlinmap_asb
end interface end interface
interface psb_is_owned
module procedure psb_is_owned
end interface
interface psb_is_local
module procedure psb_is_local
end interface
interface psb_owned_index
module procedure psb_owned_index, psb_owned_index_v
end interface
interface psb_local_index
module procedure psb_local_index, psb_local_index_v
end interface
contains contains
function psb_is_owned(idx,desc)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx
type(psb_desc_type), intent(in) :: desc
logical :: psb_is_owned
logical :: res
integer :: info
call psb_owned_index(res,idx,desc,info)
if (info /= 0) res=.false.
psb_is_owned = res
end function psb_is_owned
function psb_is_local(idx,desc)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx
type(psb_desc_type), intent(in) :: desc
logical :: psb_is_local
logical :: res
integer :: info
call psb_local_index(res,idx,desc,info)
if (info /= 0) res=.false.
psb_is_local = res
end function psb_is_local
subroutine psb_owned_index(res,idx,desc,info)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx
type(psb_desc_type), intent(in) :: desc
logical, intent(out) :: res
integer, intent(out) :: info
integer :: lx
call psb_glob_to_loc(idx,lx,desc,info,iact='I',owned=.true.)
res = (lx>0)
end subroutine psb_owned_index
subroutine psb_owned_index_v(res,idx,desc,info)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx(:)
type(psb_desc_type), intent(in) :: desc
logical, intent(out) :: res(:)
integer, intent(out) :: info
integer, allocatable :: lx(:)
allocate(lx(size(idx)),stat=info)
res=.false.
if (info /= 0) return
call psb_glob_to_loc(idx,lx,desc,info,iact='I',owned=.true.)
res = (lx>0)
end subroutine psb_owned_index_v
subroutine psb_local_index(res,idx,desc,info)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx
type(psb_desc_type), intent(in) :: desc
logical, intent(out) :: res
integer, intent(out) :: info
integer :: lx
call psb_glob_to_loc(idx,lx,desc,info,iact='I',owned=.false.)
res = (lx>0)
end subroutine psb_local_index
subroutine psb_local_index_v(res,idx,desc,info)
use psb_descriptor_type
implicit none
integer, intent(in) :: idx(:)
type(psb_desc_type), intent(in) :: desc
logical, intent(out) :: res(:)
integer, intent(out) :: info
integer, allocatable :: lx(:)
allocate(lx(size(idx)),stat=info)
res=.false.
if (info /= 0) return
call psb_glob_to_loc(idx,lx,desc,info,iact='I',owned=.false.)
res = (lx>0)
end subroutine psb_local_index_v
subroutine psb_get_boundary(bndel,desc,info) subroutine psb_get_boundary(bndel,desc,info)
use psb_descriptor_type use psb_descriptor_type
use psi_mod use psi_mod

@ -494,6 +494,10 @@ contains
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(:) integer, intent(in), optional :: ng,gtl(:)
integer :: i,ir,ic integer :: i,ir,ic
character(len=20) :: name, ch_err
name='psb_inner_upd'
if (present(gtl)) then if (present(gtl)) then
if (.not.present(ng)) then if (.not.present(ng)) then
@ -504,6 +508,7 @@ contains
do i=1, nz do i=1, nz
nza = nza + 1 nza = nza + 1
if (nza>maxsz) then if (nza>maxsz) then
call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /))
info = -71 info = -71
return return
endif endif

@ -62,11 +62,11 @@ contains
type(psb_dspmat_type), intent(inout) :: a type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
real(psb_dpk_), intent(in) :: val(*) real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
info = 0 info = 0
@ -107,11 +107,11 @@ contains
type(psb_zspmat_type), intent(inout) :: a type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
complex(psb_dpk_), intent(in) :: val(*) complex(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
info = 0 info = 0
@ -170,11 +170,11 @@ contains
type(psb_dspmat_type), intent(inout) :: a type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
real(psb_dpk_), intent(in) :: val(*) real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
character(len=20) :: name='d_csr_srch_upd' character(len=20) :: name='d_csr_srch_upd'
@ -211,51 +211,20 @@ contains
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc=i2-i1 nc=i2-i1
call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (.true.) then if (ip>0) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) a%aspk(i1+ip-1) = val(i)
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
else else
!!$ if (debug_level >= psb_debug_serial_) &
ip = -1 & write(debug_unit,*) trim(name),&
lb = i1 & ': Was searching ',ic,' in: ',i1,i2,&
ub = i2-1 & ' : ',a%ia1(i1:i2-1)
do info = i
if (lb > ub) exit return
m = (lb+ub)/2
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if end if
else else
if (debug_level >= psb_debug_serial_) & if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),& & write(debug_unit,*) trim(name),&
& ': Discarding row that does not belong to us.' & ': Discarding row that does not belong to us.'
@ -277,7 +246,7 @@ contains
i1 = a%ia2(ir) i1 = a%ia2(ir)
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc = i2-i1 nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else else
@ -322,48 +291,18 @@ contains
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc=i2-i1 nc=i2-i1
call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (.true.) then if (ip>0) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) a%aspk(i1+ip-1) = val(i)
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
else else
ip = -1 if (debug_level >= psb_debug_serial_) &
lb = i1 & write(debug_unit,*) trim(name),&
ub = i2-1 & ': Was searching ',ic,' in: ',i1,i2,&
do & ' : ',a%ia1(i1:i2-1)
if (lb > ub) exit info = i
m = (lb+ub)/2 return
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if end if
else else
if (debug_level >= psb_debug_serial_) & if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),& & write(debug_unit,*) trim(name),&
@ -383,7 +322,7 @@ contains
i1 = a%ia2(ir) i1 = a%ia2(ir)
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc = i2-i1 nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else else
@ -419,11 +358,11 @@ contains
type(psb_dspmat_type), intent(inout) :: a type(psb_dspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
real(psb_dpk_), intent(in) :: val(*) real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer :: i,ir,ic, ilr, ilc, ip, & integer :: i,ir,ic, ilr, ilc, ip, &
& i1,i2,nc,nnz,dupl & i1,i2,nc,nnz,dupl
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
@ -644,11 +583,11 @@ contains
type(psb_dspmat_type), intent(inout), target :: a type(psb_dspmat_type), intent(inout), target :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
real(psb_dpk_), intent(in) :: val(*) real(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer, pointer :: ia1(:), ia2(:), ia3(:),& integer, pointer :: ia1(:), ia2(:), ia3(:),&
& ja_(:), ka_(:) & ja_(:), ka_(:)
@ -882,11 +821,11 @@ contains
type(psb_zspmat_type), intent(inout) :: a type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
complex(psb_dpk_), intent(in) :: val(*) complex(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer :: i,ir,ic, ilr, ilc, ip, & integer :: i,ir,ic, ilr, ilc, ip, &
& i1,i2,nc,lb,ub,m,dupl & i1,i2,nc,lb,ub,m,dupl
@ -924,50 +863,18 @@ contains
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc=i2-i1 nc=i2-i1
call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (.true.) then if (ip>0) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) a%aspk(i1+ip-1) = val(i)
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
else else
!!$ if (debug_level >= psb_debug_serial_) &
ip = -1 & write(debug_unit,*) trim(name),&
lb = i1 & ': Was searching ',ic,' in: ',i1,i2,&
ub = i2-1 & ' : ',a%ia1(i1:i2-1)
do info = i
if (lb > ub) exit return
m = (lb+ub)/2
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if end if
else else
if (debug_level >= psb_debug_serial_) & if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),& & write(debug_unit,*) trim(name),&
@ -1035,48 +942,18 @@ contains
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc=i2-i1 nc=i2-i1
call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (.true.) then if (ip>0) then
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) a%aspk(i1+ip-1) = val(i)
if (ip>0) then
a%aspk(i1+ip-1) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
else else
ip = -1 if (debug_level >= psb_debug_serial_) &
lb = i1 & write(debug_unit,*) trim(name),&
ub = i2-1 & ': Was searching ',ic,' in: ',i1,i2,&
do & ' : ',a%ia1(i1:i2-1)
if (lb > ub) exit info = i
m = (lb+ub)/2 return
if (ic == a%ia1(m)) then
ip = m
lb = ub + 1
else if (ic < a%ia1(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
if (ip>0) then
a%aspk(ip) = val(i)
else
if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),&
& ': Was searching ',ic,' in: ',i1,i2,&
& ' : ',a%ia1(i1:i2-1)
info = i
return
end if
end if end if
else else
if (debug_level >= psb_debug_serial_) & if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),& & write(debug_unit,*) trim(name),&
@ -1096,7 +973,7 @@ contains
i1 = a%ia2(ir) i1 = a%ia2(ir)
i2 = a%ia2(ir+1) i2 = a%ia2(ir+1)
nc = i2-i1 nc = i2-i1
call issrch(ip,ic,nc,a%ia1(i1:i2-1)) call ibsrch(ip,ic,nc,a%ia1(i1:i2-1))
if (ip>0) then if (ip>0) then
a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i) a%aspk(i1+ip-1) = a%aspk(i1+ip-1) + val(i)
else else
@ -1132,11 +1009,11 @@ contains
type(psb_zspmat_type), intent(inout) :: a type(psb_zspmat_type), intent(inout) :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
complex(psb_dpk_), intent(in) :: val(*) complex(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer :: i,ir,ic, ilr, ilc, ip, & integer :: i,ir,ic, ilr, ilc, ip, &
& i1,i2,nc,nnz,dupl & i1,i2,nc,nnz,dupl
integer :: debug_level, debug_unit integer :: debug_level, debug_unit
@ -1360,11 +1237,11 @@ contains
type(psb_zspmat_type), intent(inout), target :: a type(psb_zspmat_type), intent(inout), target :: a
integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl integer, intent(in) :: nz, imin,imax,jmin,jmax,nzl
integer, intent(in) :: ia(*),ja(*) integer, intent(in) :: ia(:),ja(:)
integer, intent(inout) :: nza integer, intent(inout) :: nza
complex(psb_dpk_), intent(in) :: val(*) complex(psb_dpk_), intent(in) :: val(:)
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(*) integer, intent(in), optional :: ng,gtl(:)
integer, pointer :: ia1(:), ia2(:), ia3(:),& integer, pointer :: ia1(:), ia2(:), ia3(:),&
& ja_(:), ka_(:) & ja_(:), ka_(:)

@ -29,9 +29,30 @@
!!$ POSSIBILITY OF SUCH DAMAGE. !!$ POSSIBILITY OF SUCH DAMAGE.
!!$ !!$
!!$ !!$
! File: psbzcoins.f90 ! File: psb_zcoins.f90
! Subroutine: ! Subroutine: psb_zcoins
! Parameters: ! Takes a cloud of coefficients and inserts them into a sparse matrix.
! This subroutine is the serial, inner counterpart to the outer, user-level
! psb_spins.
!
! Arguments:
!
! nz - integer, input The number of points to insert.
! ia(:) - integer, input The row indices of the coefficients.
! ja(:) - integer, input The column indices of the coefficients.
! val(:) - complex, input The values of the coefficients to be inserted.
! a - type(psb_zspmat_type), inout The sparse destination matrix.
! imin - integer, input The minimum valid row index
! imax - integer, input The maximum valid row index
! jmin - integer, input The minimum valid col index
! jmax - integer, input The maximum valid col index
! info - integer, output Return code.
! gtl(:) - integer, input,optional An index mapping to be applied
! default: identity
! rebuild - logical, input, optional Rebuild in case of update
! finding a new index. Default: false.
! Not fully tested.
!
subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild) subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
use psb_spmat_type use psb_spmat_type
@ -54,7 +75,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
character(len=5) :: ufida character(len=5) :: ufida
integer :: ng, nza, isza,spstate, & integer :: ng, nza, isza,spstate, &
& ip1, nzl, err_act, int_err(5), iupd, irst & ip1, nzl, err_act, int_err(5), iupd, irst
logical, parameter :: debug=.false. integer :: debug_level, debug_unit
logical :: rebuild_ logical :: rebuild_
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
type(psb_zspmat_type) :: tmp type(psb_zspmat_type) :: tmp
@ -62,6 +83,8 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
name='psb_zcoins' name='psb_zcoins'
info = 0 info = 0
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
info = 0 info = 0
if (nz <= 0) then if (nz <= 0) then
@ -140,9 +163,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999 goto 9999
endif endif
if (debug) then if (debug_level >= psb_debug_serial_) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items ' write(debug_unit,*) trim(name),': insert discarded items '
end if end if
end if end if
if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then
@ -172,9 +195,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
goto 9999 goto 9999
endif endif
if (debug) then if (debug_level >= psb_debug_serial_) then
if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: update discarded items ' write(debug_unit,*) trim(name),': update discarded items '
end if end if
end if end if
@ -185,7 +208,8 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
a%ia2(ip1+psb_nnz_) = nza a%ia2(ip1+psb_nnz_) = nza
end select end select
if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza if (debug_level >= psb_debug_serial_) &
& write(debug_unit,*) trim(name),': (UPD) : NZA:',nza
case (psb_upd_srch_) case (psb_upd_srch_)
@ -194,8 +218,10 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
if (info > 0) then if (info > 0) then
if (rebuild_) then if (rebuild_) then
if (debug) write(0,*)& if (debug_level >= psb_debug_serial_) &
& 'COINS: Going through rebuild_ fingers crossed!' & write(debug_unit,*)&
& trim(name),&
& ': Going through rebuild_ fingers crossed!'
irst = info irst = info
call psb_nullify_sp(tmp) call psb_nullify_sp(tmp)
call psb_spcnv(a,tmp,info,afmt='coo') call psb_spcnv(a,tmp,info,afmt='coo')
@ -206,9 +232,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
goto 9999 goto 9999
endif endif
call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info)
if (debug) then if (debug_level >= psb_debug_serial_) &
write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst & write(debug_unit,*) trim(name),&
endif & ': Rebuild size',tmp%infoa(psb_nnz_) ,irst
call psb_sp_transfer(tmp,a,info) call psb_sp_transfer(tmp,a,info)
if(info /= izero) then if(info /= izero) then
info=4010 info=4010
@ -225,8 +251,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
goto 9999 goto 9999
endif endif
if (debug) write(0,*)& if (debug_level >= psb_debug_serial_) write(debug_unit,*)&
& 'COINS: Reinserting',a%fida,nza,isza,irst,nz & trim(name),': Reinserting',a%fida,nza,isza,irst,nz
if ((nza+nz)>isza) then if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info /= izero) then if(info /= izero) then
@ -245,7 +272,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
ch_err='psb_inner_ins' ch_err='psb_inner_ins'
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
endif endif
call psb_sp_setifld(nza,psb_del_bnd_,a,info) call psb_sp_setifld(nza,psb_del_bnd_,a,info)
call psb_sp_setifld(nza,psb_nnz_,a,info) call psb_sp_setifld(nza,psb_nnz_,a,info)
end if end if
@ -321,14 +348,13 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
goto 9999 goto 9999
endif endif
if (debug) then if (debug_level >= psb_debug_serial_) then
if ((nza - a%infoa(psb_nnz_)) /= nz) then if ((nza - a%infoa(psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: insert discarded items ' write(debug_unit,*) trim(name),': insert discarded items '
end if end if
end if end if
if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then if ((nza - psb_sp_getifld(psb_nnz_,a,info)) /= nz) then
call psb_sp_setifld(nza,psb_del_bnd_,a,info) call psb_sp_setifld(nza,psb_del_bnd_,a,info)
!!$ write(0,*) 'Settind del_bnd_ 2: ',nza
endif endif
call psb_sp_setifld(nza,psb_nnz_,a,info) call psb_sp_setifld(nza,psb_nnz_,a,info)
@ -350,14 +376,15 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
goto 9999 goto 9999
endif endif
if (debug) then if (debug_level >= psb_debug_serial_) then
if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then if ((nza - a%ia2(ip1+psb_nnz_)) /= nz) then
write(0,*) 'PSB_COINS: update discarded items ' write(debug_unit,*) trim(name),': update discarded items '
end if end if
end if end if
a%ia2(ip1+psb_nnz_) = nza a%ia2(ip1+psb_nnz_) = nza
if (debug) write(0,*) 'From COINS(UPD) : NZA:',nza if (debug_level >= psb_debug_serial_)&
& write(debug_unit,*) trim(name),':(UPD) : NZA:',nza
case (psb_upd_srch_) case (psb_upd_srch_)
@ -366,15 +393,17 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
if (info > 0) then if (info > 0) then
if (rebuild_) then if (rebuild_) then
if (debug) write(0,*)& if (debug_level >= psb_debug_serial_) &
& 'COINS: Going through rebuild_ fingers crossed!' & write(debug_unit,*)&
& trim(name),&
& ': Going through rebuild_ fingers crossed!'
irst = info irst = info
call psb_nullify_sp(tmp) call psb_nullify_sp(tmp)
call psb_spcnv(a,tmp,info,afmt='coo') call psb_spcnv(a,tmp,info,afmt='coo')
call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info) call psb_sp_setifld(psb_spmat_bld_,psb_state_,tmp,info)
if (debug) then if (debug_level >= psb_debug_serial_) &
write(0,*) 'COINS Rebuild: size',tmp%infoa(psb_nnz_) ,irst & write(debug_unit,*) trim(name),&
endif & ': Rebuild size',tmp%infoa(psb_nnz_) ,irst
call psb_sp_transfer(tmp,a,info) call psb_sp_transfer(tmp,a,info)
call psb_sp_info(psb_nztotreq_,a,nza,info) call psb_sp_info(psb_nztotreq_,a,nza,info)
call psb_sp_info(psb_nzsizereq_,a,isza,info) call psb_sp_info(psb_nzsizereq_,a,isza,info)
@ -385,8 +414,9 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
goto 9999 goto 9999
endif endif
if (debug) write(0,*)& if (debug_level >= psb_debug_serial_) write(debug_unit,*)&
& 'COINS: Reinserting',a%fida,nza,isza & trim(name),': Reinserting',a%fida,nza,isza
if ((nza+nz)>isza) then if ((nza+nz)>isza) then
call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info) call psb_sp_reall(a,max(nza+nz,int(1.5*isza)),info)
if(info /= izero) then if(info /= izero) then
@ -405,7 +435,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
ch_err='psb_inner_ins' ch_err='psb_inner_ins'
call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/)) call psb_errpush(4013,name,a_err=ch_err,i_err=(/info,0,0,0,0/))
endif endif
call psb_sp_setifld(nza,psb_del_bnd_,a,info) call psb_sp_setifld(nza,psb_del_bnd_,a,info)
call psb_sp_setifld(nza,psb_nnz_,a,info) call psb_sp_setifld(nza,psb_nnz_,a,info)
end if end if
@ -438,6 +468,7 @@ subroutine psb_zcoins(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl,rebuild)
endif endif
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -463,6 +494,10 @@ contains
integer, intent(out) :: info integer, intent(out) :: info
integer, intent(in), optional :: ng,gtl(:) integer, intent(in), optional :: ng,gtl(:)
integer :: i,ir,ic integer :: i,ir,ic
character(len=20) :: name, ch_err
name='psb_inner_upd'
if (present(gtl)) then if (present(gtl)) then
if (.not.present(ng)) then if (.not.present(ng)) then
@ -473,7 +508,7 @@ contains
do i=1, nz do i=1, nz
nza = nza + 1 nza = nza + 1
if (nza>maxsz) then if (nza>maxsz) then
write(0,*) 'Out of bounds in INNER_UPD ',nza,maxsz call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /))
info = -71 info = -71
return return
endif endif

@ -1058,11 +1058,6 @@ An integer value; 0 means no error has been detected.
\item The default \verb|I|gnore means that the negative output is the \item The default \verb|I|gnore means that the negative output is the
only action taken on an out-of-range input. only action taken on an out-of-range input.
\end{enumerate} \end{enumerate}
%
%% psb_loc_to_glob %%
%
\subroutine{psb\_loc\_to\_glob}{Local to global indices conversion} \subroutine{psb\_loc\_to\_glob}{Local to global indices conversion}
\syntax{call psb\_loc\_to\_glob}{x, y, desc\_a, info, iact} \syntax{call psb\_loc\_to\_glob}{x, y, desc\_a, info, iact}
@ -1113,6 +1108,180 @@ An integer value; 0 means no error has been detected.
%
%% psb_loc_to_glob %%
%
\subroutine{psb\_is\_owned}{}
\syntax{call psb\_is\_owned}{x, desc\_a}
\begin{description}
\item[Type:] Asynchronous.
\item[\bf On Entry]
\item[x] Integer index.\\
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf in}.\\
Specified as: a scalar integer.\\
\item[desc\_a] the communication descriptor.\\
Scope:{\bf local}.\\
Type:{\bf required}.\\
Intent: {\bf in}.\\
Specified as: a structured data of type \descdata.
\end{description}
\begin{description}
\item[\bf On Return]
\item[Function value] A logical mask which is true if
$x$ is owned by the current process
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf out}.\\
\end{description}
\section*{Notes}
\begin{enumerate}
\item This routine returns a \verb|.true.| value for an index
that is strictly owned by the current process, excluding the halo
indices
\end{enumerate}
\subroutine{psb\_owned\_index}{}
\syntax{call psb\_owned\_index}{y, x, desc\_a, info}
\begin{description}
\item[Type:] Asynchronous.
\item[\bf On Entry]
\item[x] Integer indices.\\
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf in, inout}.\\
Specified as: a scalar or a rank one integer array.\\
\item[desc\_a] the communication descriptor.\\
Scope:{\bf local}.\\
Type:{\bf required}.\\
Intent: {\bf in}.\\
Specified as: a structured data of type \descdata.
\item[iact] specifies action to be taken in case of range errors.
Scope: {\bf global} \\
Type: {\bf optional}\\
Intent: {\bf in}.\\
Specified as: a character variable \verb|I|gnore, \verb|W|arning or
\verb|A|bort, default \verb|I|gnore.
\end{description}
\begin{description}
\item[\bf On Return]
\item[y] A logical mask which is true for all corresponding entries of
$x$ that are owned by the current process
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf out}.\\
Specified as: a scalar or rank one logical array.
\item[info] Error code.\\
Scope: {\bf local} \\
Type: {\bf required} \\
Intent: {\bf out}.\\
An integer value; 0 means no error has been detected.
\end{description}
\section*{Notes}
\begin{enumerate}
\item This routine returns a \verb|.true.| value for those indices
that are strictly owned by the current process, excluding the halo
indices
\end{enumerate}
\subroutine{psb\_is\_local}{}
\syntax{call psb\_is\_local}{x, desc\_a}
\begin{description}
\item[Type:] Asynchronous.
\item[\bf On Entry]
\item[x] Integer index.\\
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf in}.\\
Specified as: a scalar integer.\\
\item[desc\_a] the communication descriptor.\\
Scope:{\bf local}.\\
Type:{\bf required}.\\
Intent: {\bf in}.\\
Specified as: a structured data of type \descdata.
\end{description}
\begin{description}
\item[\bf On Return]
\item[Function value] A logical mask which is true if
$x$ is local to the current process
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf out}.\\
\end{description}
\section*{Notes}
\begin{enumerate}
\item This routine returns a \verb|.true.| value for an index
that is local to the current process, including the halo
indices
\end{enumerate}
\subroutine{psb\_local\_index}{}
\syntax{call psb\_local\_index}{y, x, desc\_a, info}
\begin{description}
\item[Type:] Asynchronous.
\item[\bf On Entry]
\item[x] Integer indices.\\
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf in, inout}.\\
Specified as: a scalar or a rank one integer array.\\
\item[desc\_a] the communication descriptor.\\
Scope:{\bf local}.\\
Type:{\bf required}.\\
Intent: {\bf in}.\\
Specified as: a structured data of type \descdata.
\item[iact] specifies action to be taken in case of range errors.
Scope: {\bf global} \\
Type: {\bf optional}\\
Intent: {\bf in}.\\
Specified as: a character variable \verb|I|gnore, \verb|W|arning or
\verb|A|bort, default \verb|I|gnore.
\end{description}
\begin{description}
\item[\bf On Return]
\item[y] A logical mask which is true for all corresponding entries of
$x$ that are local to the current process
Scope: {\bf local} \\
Type: {\bf required}\\
Intent: {\bf out}.\\
Specified as: a scalar or rank one logical array.
\item[info] Error code.\\
Scope: {\bf local} \\
Type: {\bf required} \\
Intent: {\bf out}.\\
An integer value; 0 means no error has been detected.
\end{description}
\section*{Notes}
\begin{enumerate}
\item This routine returns a \verb|.true.| value for those indices
that are local to the current process, including the halo
indices.
\end{enumerate}
% %
%% psb_ins %% %% psb_ins %%

File diff suppressed because one or more lines are too long

@ -19,6 +19,10 @@ ppde: ppde.o
$(F90LINK) ppde.o -o ppde $(PSBLAS_LIB) $(LDLIBS) $(F90LINK) ppde.o -o ppde $(PSBLAS_LIB) $(LDLIBS)
/bin/mv ppde $(EXEDIR) /bin/mv ppde $(EXEDIR)
tpde: tpde.o
$(F90LINK) tpde.o -pg -o tpde $(PSBLAS_LIB) $(LDLIBS)
/bin/mv tpde $(EXEDIR)
.f90.o: .f90.o:
$(MPF90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $< $(MPF90) $(F90COPT) $(FINCLUDES) $(FDEFINES) -c $<

@ -361,14 +361,12 @@ contains
end interface ! local variables end interface ! local variables
type(psb_dspmat_type) :: a type(psb_dspmat_type) :: a
real(psb_dpk_) :: zt(nbmax),glob_x,glob_y,glob_z real(psb_dpk_) :: zt(nbmax),glob_x,glob_y,glob_z
integer :: m,n,nnz,glob_row integer :: m,n,nnz,glob_row,loc_row
integer :: x,y,z,ia,indx_owner integer :: x,y,z,ia,indx_owner
integer :: np, iam integer :: np, iam
integer :: element integer :: element
integer :: nv, inv
integer, allocatable :: irow(:),icol(:) integer, allocatable :: irow(:),icol(:)
real(psb_dpk_), allocatable :: val(:) real(psb_dpk_), allocatable :: val(:)
integer, allocatable :: prv(:)
! deltah dimension of each grid cell ! deltah dimension of each grid cell
! deltat discretization time ! deltat discretization time
real(psb_dpk_) :: deltah real(psb_dpk_) :: deltah
@ -414,7 +412,7 @@ contains
! a bunch of rows per call. ! a bunch of rows per call.
! !
allocate(val(20*nbmax),irow(20*nbmax),& allocate(val(20*nbmax),irow(20*nbmax),&
&icol(20*nbmax),prv(np),stat=info) &icol(20*nbmax),stat=info)
if (info /= 0 ) then if (info /= 0 ) then
info=4000 info=4000
call psb_errpush(info,name) call psb_errpush(info,name)
@ -430,138 +428,135 @@ contains
! icol(1)=1 ! icol(1)=1
do glob_row = 1, n do glob_row = 1, n
call parts(glob_row,n,np,prv,nv) ! Figure out which rows are local to the current process:
do inv = 1, nv if (psb_is_owned(glob_row,desc_a)) then
indx_owner = prv(inv) ! local matrix pointer
if (indx_owner == iam) then element=1
! local matrix pointer ! compute gridpoint coordinates
element=1 if (mod(glob_row,(idim*idim)) == 0) then
! compute gridpoint coordinates x = glob_row/(idim*idim)
if (mod(glob_row,(idim*idim)) == 0) then else
x = glob_row/(idim*idim) x = glob_row/(idim*idim)+1
else endif
x = glob_row/(idim*idim)+1 if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then
endif y = (glob_row-(x-1)*idim*idim)/idim
if (mod((glob_row-(x-1)*idim*idim),idim) == 0) then else
y = (glob_row-(x-1)*idim*idim)/idim y = (glob_row-(x-1)*idim*idim)/idim+1
else endif
y = (glob_row-(x-1)*idim*idim)/idim+1 z = glob_row-(x-1)*idim*idim-(y-1)*idim
endif ! glob_x, glob_y, glob_x coordinates
z = glob_row-(x-1)*idim*idim-(y-1)*idim glob_x=x*deltah
! glob_x, glob_y, glob_x coordinates glob_y=y*deltah
glob_x=x*deltah glob_z=z*deltah
glob_y=y*deltah
glob_z=z*deltah ! check on boundary points
zt(1) = 0.d0
! check on boundary points ! internal point: build discretization
zt(1) = 0.d0 !
! internal point: build discretization ! term depending on (x-1,y,z)
! !
! term depending on (x-1,y,z) if (x==1) then
! val(element)=-b1(glob_x,glob_y,glob_z)&
if (x==1) then & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element))
& deltah) else
zt(1) = exp(-glob_y**2-glob_z**2)*(-val(element)) val(element)=-b1(glob_x,glob_y,glob_z)&
else & -a1(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z)& val(element) = val(element)/(deltah*&
& -a1(glob_x,glob_y,glob_z) & deltah)
val(element) = val(element)/(deltah*& icol(element)=(x-2)*idim*idim+(y-1)*idim+(z)
& deltah) element=element+1
icol(element)=(x-2)*idim*idim+(y-1)*idim+(z) endif
element=element+1 ! term depending on (x,y-1,z)
endif if (y==1) then
! term depending on (x,y-1,z) val(element)=-b2(glob_x,glob_y,glob_z)&
if (y==1) then & -a2(glob_x,glob_y,glob_z)
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b2(glob_x,glob_y,glob_z)&
& -a2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
element=element+1
endif
! term depending on (x,y,z-1)
if (z==1) then
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b3(glob_x,glob_y,glob_z)&
& -a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
element=element+1
endif
! term depending on (x,y,z)
val(element)=2*b1(glob_x,glob_y,glob_z)&
& +2*b2(glob_x,glob_y,glob_z)&
& +2*b3(glob_x,glob_y,glob_z)&
& +a1(glob_x,glob_y,glob_z)&
& +a2(glob_x,glob_y,glob_z)&
& +a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
element=element+1 else
! term depending on (x,y,z+1) val(element)=-b2(glob_x,glob_y,glob_z)&
if (z==idim) then & -a2(glob_x,glob_y,glob_z)
val(element)=-b1(glob_x,glob_y,glob_z) val(element) = val(element)/(deltah*&
val(element) = val(element)/(deltah*& & deltah)
& deltah) icol(element)=(x-1)*idim*idim+(y-2)*idim+(z)
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) element=element+1
else endif
val(element)=-b1(glob_x,glob_y,glob_z) ! term depending on (x,y,z-1)
val(element) = val(element)/(deltah*& if (z==1) then
& deltah) val(element)=-b3(glob_x,glob_y,glob_z)&
icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1) & -a3(glob_x,glob_y,glob_z)
element=element+1 val(element) = val(element)/(deltah*&
endif & deltah)
! term depending on (x,y+1,z) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
if (y==idim) then else
val(element)=-b2(glob_x,glob_y,glob_z) val(element)=-b3(glob_x,glob_y,glob_z)&
val(element) = val(element)/(deltah*& & -a3(glob_x,glob_y,glob_z)
& deltah) val(element) = val(element)/(deltah*&
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element)) & deltah)
else icol(element)=(x-1)*idim*idim+(y-1)*idim+(z-1)
val(element)=-b2(glob_x,glob_y,glob_z) element=element+1
val(element) = val(element)/(deltah*& endif
& deltah) ! term depending on (x,y,z)
icol(element)=(x-1)*idim*idim+(y)*idim+(z) val(element)=2*b1(glob_x,glob_y,glob_z)&
element=element+1 & +2*b2(glob_x,glob_y,glob_z)&
endif & +2*b3(glob_x,glob_y,glob_z)&
! term depending on (x+1,y,z) & +a1(glob_x,glob_y,glob_z)&
if (x<idim) then & +a2(glob_x,glob_y,glob_z)&
val(element)=-b3(glob_x,glob_y,glob_z) & +a3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*& val(element) = val(element)/(deltah*&
& deltah) & deltah)
icol(element)=(x)*idim*idim+(y-1)*idim+(z) icol(element)=(x-1)*idim*idim+(y-1)*idim+(z)
element=element+1 element=element+1
endif ! term depending on (x,y,z+1)
irow(1:element-1)=glob_row if (z==idim) then
ia=glob_row val(element)=-b1(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
t3 = psb_wtime() & deltah)
call psb_spins(element-1,irow,icol,val,a,desc_a,info) zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
if(info /= 0) exit else
tins = tins + (psb_wtime()-t3) val(element)=-b1(glob_x,glob_y,glob_z)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info) val(element) = val(element)/(deltah*&
if(info /= 0) exit & deltah)
zt(1)=0.d0 icol(element)=(x-1)*idim*idim+(y-1)*idim+(z+1)
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info) element=element+1
if(info /= 0) exit endif
end if ! term depending on (x,y+1,z)
end do if (y==idim) then
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
zt(1) = exp(-glob_y**2-glob_z**2)*exp(-glob_x)*(-val(element))
else
val(element)=-b2(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x-1)*idim*idim+(y)*idim+(z)
element=element+1
endif
! term depending on (x+1,y,z)
if (x<idim) then
val(element)=-b3(glob_x,glob_y,glob_z)
val(element) = val(element)/(deltah*&
& deltah)
icol(element)=(x)*idim*idim+(y-1)*idim+(z)
element=element+1
endif
irow(1:element-1)=glob_row
ia=glob_row
t3 = psb_wtime()
call psb_spins(element-1,irow,icol,val,a,desc_a,info)
if(info /= 0) exit
tins = tins + (psb_wtime()-t3)
call psb_geins(1,(/ia/),zt(1:1),b,desc_a,info)
if(info /= 0) exit
zt(1)=0.d0
call psb_geins(1,(/ia/),zt(1:1),xv,desc_a,info)
if(info /= 0) exit
end if
end do end do
call psb_barrier(ictxt) call psb_barrier(ictxt)

Loading…
Cancel
Save