psblas3-testmv:

base/modules/Makefile
 base/modules/psb_d_parmat_mod.f90
 base/modules/psb_d_psblas_mod.f90
 base/modules/psb_d_tools_mod.f90
 base/psblas/psb_dspmm.f90
 base/tools/psb_dspalloc.f90
 base/tools/psb_dspasb.f90
 base/tools/psb_dspfree.f90
 krylov/psb_dcgstab.f90
 krylov/psb_dkrylov.f90
 krylov/psb_krylov_mod.f90
 test/kernel/Makefile
 test/kernel/pgenpmv.f90
 test/kernel/runs/spmv.inp
 test/pargen/runs/ppde.inp
 util/psb_d_genpde_impl.f90
 util/psb_d_genpde_mod.f90

Try to overlap comm and comp.
psblas-testmv
Salvatore Filippone 11 years ago
parent 32f5f86e9e
commit 39ed3311e7

@ -27,6 +27,7 @@ UTIL_MODS = psb_string_mod.o psb_desc_const_mod.o psb_indx_map_mod.o\
psb_ip_reord_mod.o\
psb_check_mod.o psb_hash_mod.o\
psb_base_mat_mod.o psb_mat_mod.o\
psb_d_parmat_mod.o \
psb_s_base_mat_mod.o psb_s_csr_mat_mod.o psb_s_csc_mat_mod.o psb_s_mat_mod.o \
psb_d_base_mat_mod.o psb_d_csr_mat_mod.o psb_d_csc_mat_mod.o psb_d_mat_mod.o \
psb_c_base_mat_mod.o psb_c_csr_mat_mod.o psb_c_csc_mat_mod.o psb_c_mat_mod.o \
@ -134,6 +135,8 @@ psb_d_comm_mod.o: psb_d_vect_mod.o psb_desc_mod.o psb_mat_mod.o
psb_c_comm_mod.o: psb_c_vect_mod.o psb_desc_mod.o psb_mat_mod.o
psb_z_comm_mod.o: psb_z_vect_mod.o psb_desc_mod.o psb_mat_mod.o
psb_d_parmat_mod.o: psb_d_mat_mod.o psb_d_vect_mod.o psb_error_mod.o psb_desc_mod.o
psb_base_mod.o: $(MODULES)
psi_penv_mod.o: psi_penv_mod.F90 $(BASIC_MODS)

File diff suppressed because it is too large Load Diff

@ -294,7 +294,7 @@ module psb_d_psblas_mod
& desc_a, info, trans, work,doswap)
import :: psb_desc_type, psb_dpk_, psb_ipk_, &
& psb_d_vect_type, psb_dspmat_type
type(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(in) :: a
real(psb_dpk_), intent(inout), target :: x(:)
real(psb_dpk_), intent(inout), target :: y(:)
real(psb_dpk_), intent(in) :: alpha, beta
@ -308,7 +308,7 @@ module psb_d_psblas_mod
& desc_a, info, trans, work,doswap)
import :: psb_desc_type, psb_dpk_, psb_ipk_, &
& psb_d_vect_type, psb_dspmat_type
type(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(in) :: a
type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y
real(psb_dpk_), intent(in) :: alpha, beta

@ -241,7 +241,7 @@ Module psb_d_tools_mod
& psb_d_base_vect_type, psb_d_vect_type, &
& psb_dspmat_type, psb_d_base_sparse_mat
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), intent(inout) :: a
class(psb_dspmat_type), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: nnz
end subroutine psb_dspalloc
@ -252,7 +252,7 @@ Module psb_d_tools_mod
import :: psb_desc_type, psb_dpk_, psb_ipk_, &
& psb_d_base_vect_type, psb_d_vect_type, &
& psb_dspmat_type, psb_d_base_sparse_mat
type(psb_dspmat_type), intent (inout) :: a
class(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd
@ -267,7 +267,7 @@ Module psb_d_tools_mod
& psb_d_base_vect_type, psb_d_vect_type, &
& psb_dspmat_type, psb_d_base_sparse_mat
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), intent(inout) ::a
class(psb_dspmat_type), intent(inout) ::a
integer(psb_ipk_), intent(out) :: info
end subroutine psb_dspfree
end interface

@ -424,8 +424,8 @@ subroutine psb_dspmv(alpha,a,x,beta,y,desc_a,info,&
real(psb_dpk_), intent(in) :: alpha, beta
real(psb_dpk_), intent(inout), target :: x(:)
real(psb_dpk_), intent(inout), target :: y(:)
type(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
class(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional, target, intent(inout) :: work(:)
character, intent(in), optional :: trans
@ -685,7 +685,7 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
real(psb_dpk_), intent(in) :: alpha, beta
type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y
type(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
real(psb_dpk_), optional, target, intent(inout) :: work(:)
@ -886,3 +886,203 @@ subroutine psb_dspmv_vect(alpha,a,x,beta,y,desc_a,info,&
end if
return
end subroutine psb_dspmv_vect
subroutine psb_d_p_csmv_vect(alpha,a,x,beta,y,info,trans)
use psb_base_mod
use psb_d_parmat_mod, psb_protect_name => psb_d_p_csmv_vect
use psi_mod
implicit none
class(psb_d_parmat_type), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta
type(psb_d_vect_type), intent(inout) :: x
type(psb_d_vect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
! locals
integer(psb_ipk_) :: ictxt, np, me,&
& err_act, n, iix, jjx, ia, ja, iia, jja, ix, iy, ik, &
& m, nrow, ncol, lldx, lldy, liwork, jx, jy, iiy, jjy,&
& ib, ip, idx
integer(psb_ipk_), parameter :: nb=4
real(psb_dpk_), pointer :: iwork(:), xp(:), yp(:)
real(psb_dpk_), allocatable :: xvsave(:)
character :: trans_
character(len=20) :: name, ch_err
logical :: aliw, doswap_
integer(psb_ipk_) :: debug_level, debug_unit
name='psb_dspmv'
if (psb_errstatus_fatal()) return
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt=a%desc%get_context()
call psb_info(ictxt, me, np)
if (np == -1) then
info = psb_err_context_error_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(x%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (.not.allocated(y%v)) then
info = psb_err_invalid_vect_state_
call psb_errpush(info,name)
goto 9999
endif
if (present(trans)) then
trans_ = psb_toupper(trans)
else
trans_ = 'N'
endif
if ( (trans_ == 'N').or.(trans_ == 'T')&
& .or.(trans_ == 'C')) then
else
info = psb_err_iarg_invalid_value_
call psb_errpush(info,name)
goto 9999
end if
m = a%desc%get_global_rows()
n = a%desc%get_global_cols()
nrow = a%desc%get_local_rows()
ncol = a%desc%get_local_cols()
lldx = x%get_nrows()
lldy = y%get_nrows()
if ((info == 0).and.(lldx<ncol)) call x%reall(ncol,info)
if ((info == 0).and.(lldy<ncol)) call y%reall(ncol,info)
if (psb_errstatus_fatal()) then
info=psb_err_from_subroutine_
ch_err='reall'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
iwork => null()
! check for presence/size of a work area
liwork= 2*ncol
aliw=.true.
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='Allocate'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
endif
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' Allocated work ', info
if (trans_ == 'N') then
! Matrix is not transposed
call psi_swapdata(psb_swap_send_,&
& dzero,x%v,a%desc,iwork,info,data=psb_comm_halo_)
if (info ==0) call psb_csmm(alpha,a%a_d,x,beta,y,info)
if (info ==0) call psi_swapdata(psb_swap_recv_,&
& dzero,x%v,a%desc,iwork,info,data=psb_comm_halo_)
if (info ==0) call psb_csmm(alpha,a%a_nd,x,done,y,info)
if(info /= psb_success_) then
info = psb_err_from_subroutine_non_
call psb_errpush(info,name)
goto 9999
end if
else
! Matrix is transposed
info = psb_err_missing_override_method_
ch_err='psb_csmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
!
! Non-empty overlap, need a buffer to hold
! the entries updated with average operator.
! Why the average? because in this way they will contribute
! with a proper scale factor (1/np) to the overall product.
!
call psi_ovrl_save(x%v,xvsave,a%desc,info)
if (info == psb_success_) call psi_ovrl_upd(x%v,a%desc,psb_avg_,info)
!!! THIS SHOULD BE FIXED !!! But beta is almost never /= 0
!!$ yp(nrow+1:ncol) = dzero
! local Matrix-vector product
!!$ if (info == psb_success_) call psb_csmm(alpha,a,x,beta,y,info,trans=trans_)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' csmm ', info
if (info == psb_success_) call psi_ovrl_restore(x%v,xvsave,a%desc,info)
if (info /= psb_success_) then
info = psb_err_from_subroutine_
ch_err='psb_csmm'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
call psi_swaptran(ior(psb_swap_send_,psb_swap_recv_),&
& done,y%v,a%desc,iwork,info)
if (info == psb_success_) call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& done,y%v,a%desc,iwork,info,data=psb_comm_ovr_)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' swaptran ', info
if(info /= psb_success_) then
info = psb_err_from_subroutine_
ch_err='PSI_SwapTran'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
end if
if (aliw) deallocate(iwork,stat=info)
if (debug_level >= psb_debug_comp_) &
& write(debug_unit,*) me,' ',trim(name),' deallocat ',aliw, info
if(info /= psb_success_) then
info = psb_err_from_subroutine_
ch_err='Deallocate iwork'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
nullify(iwork)
call psb_erractionrestore(err_act)
if (debug_level >= psb_debug_comp_) then
call psb_barrier(ictxt)
write(debug_unit,*) me,' ',trim(name),' Returning '
endif
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error(ictxt)
return
end if
return
end subroutine psb_d_p_csmv_vect

@ -47,7 +47,7 @@ subroutine psb_dspalloc(a, desc_a, info, nnz)
!....parameters...
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), intent(inout) :: a
class(psb_dspmat_type), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), optional, intent(in) :: nnz

@ -55,7 +55,7 @@ subroutine psb_dspasb(a,desc_a, info, afmt, upd, dupl, mold)
!...Parameters....
type(psb_dspmat_type), intent (inout) :: a
class(psb_dspmat_type), intent (inout) :: a
type(psb_desc_type), intent(in) :: desc_a
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_),optional, intent(in) :: dupl, upd

@ -45,7 +45,7 @@ subroutine psb_dspfree(a, desc_a,info)
!....parameters...
type(psb_desc_type), intent(in) :: desc_a
type(psb_dspmat_type), intent(inout) :: a
class(psb_dspmat_type), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
!...locals....
integer(psb_ipk_) :: ictxt,err_act

@ -99,7 +99,7 @@ Subroutine psb_dcgstab_vect(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,ist
use psb_d_krylov_conv_mod
use psb_krylov_mod
implicit none
type(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(in) :: a
class(psb_dprec_type), Intent(inout) :: prec
Type(psb_desc_type), Intent(in) :: desc_a
type(psb_d_vect_type), Intent(inout) :: b

@ -86,7 +86,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
use psb_krylov_mod, psb_protect_name => psb_dkrylov_vect
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
class(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b
@ -132,7 +132,7 @@ Subroutine psb_dkrylov_vect(method,a,prec,b,x,eps,desc_a,info,&
& desc_a,info,itmax,iter,err,itrace,istop)
import :: psb_ipk_, psb_dpk_, psb_desc_type, &
& psb_dspmat_type, psb_dprec_type, psb_d_vect_type
type(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(in) :: a
type(psb_desc_type), intent(in) :: desc_a
type(psb_d_vect_type), Intent(inout) :: b
type(psb_d_vect_type), Intent(inout) :: x

@ -90,7 +90,7 @@ Module psb_krylov_mod
use psb_prec_mod, only : psb_dprec_type
character(len=*) :: method
Type(psb_dspmat_type), Intent(in) :: a
class(psb_dspmat_type), Intent(in) :: a
Type(psb_desc_type), Intent(in) :: desc_a
class(psb_dprec_type), intent(inout) :: prec
type(psb_d_vect_type), Intent(inout) :: b

@ -12,9 +12,10 @@ DBOBJS=dpsbench.o
DTOBJS=d_file_spmv.o
STOBJS=s_file_spmv.o
DPGOBJS=pdgenspmv.o
PPOBJS=pgenpmv.o
EXEDIR=./runs
all: d_file_spmv s_file_spmv pdgenspmv
all: d_file_spmv s_file_spmv pdgenspmv pgenpmv
dpsbench: $(DBOBJS)
@ -28,6 +29,9 @@ d_file_spmv: $(DTOBJS)
pdgenspmv: $(DPGOBJS)
$(F90LINK) $(LOPT) $(DPGOBJS) -o pdgenspmv $(PSBLAS_LIB) $(LDLIBS)
/bin/mv pdgenspmv $(EXEDIR)
pgenpmv: $(PPOBJS)
$(F90LINK) $(LOPT) $(PPOBJS) -o pgenpmv $(PSBLAS_LIB) $(LDLIBS)
/bin/mv pgenpmv $(EXEDIR)
s_file_spmv: $(STOBJS)

@ -0,0 +1,321 @@
!!$
!!$ Parallel Sparse BLAS version 2.3.1
!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012
!!$ Salvatore Filippone University of Rome Tor Vergata
!!$ Alfredo Buttari CNRS-IRIT, Toulouse
!!$
!!$ 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.
!!$
!!$
! File: ppde.f90
!
program pgenpmv
use psb_base_mod
use psb_d_parmat_mod
use psb_util_mod
implicit none
! input parameters
character(len=20) :: kmethd, ptype
character(len=5) :: afmt
integer(psb_ipk_) :: idim
! miscellaneous
real(psb_dpk_), parameter :: one = 1.d0
real(psb_dpk_) :: t1, t2, tprec, flops, tflops, tt1, tt2, bdwdth
! sparse matrix and preconditioner
type(psb_dspmat_type) :: a
type(psb_d_parmat_type) :: ap
! descriptor
type(psb_desc_type) :: desc_a
! dense matrices
type(psb_d_vect_type) :: xv,bv, vtst
real(psb_dpk_), allocatable :: tst(:)
! blacs parameters
integer(psb_ipk_) :: ictxt, iam, np
! solver parameters
integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst, nr
integer(psb_long_int_k_) :: amatsize, precsize, descsize, d2size, annz, nbytes
real(psb_dpk_) :: err, eps, epsi, eps2
integer(psb_ipk_), parameter :: times=10
! other variables
integer(psb_ipk_) :: info, i
character(len=20) :: name,ch_err
character(len=40) :: fname
info=psb_success_
call psb_init(ictxt)
call psb_info(ictxt,iam,np)
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ictxt)
stop
endif
if(psb_get_errstatus() /= 0) goto 9999
name='pde90'
call psb_set_errverbosity(itwo)
!
! Hello world
!
if (iam == psb_root_) then
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
end if
!
! get parameters
!
call get_parms(ictxt,afmt,idim)
!
! allocate and fill in the coefficient matrix, rhs and initial guess
!
call psb_barrier(ictxt)
t1 = psb_wtime()
call psb_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
& a1,a2,a3,b1,b2,b3,c,g,info)
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_gen_pde3d'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(psb_out_unit,'(" ")')
call ap%create(a,desc_a,info)
write(0,*) 'After ap%create: ',info,&
& ap%sizeof(),ap%a_d%get_nzeros(),ap%a_nd%get_nzeros()
call xv%set(done)
call vtst%bld(xv%get_nrows())
call psb_barrier(ictxt)
t1 = psb_wtime()
do i=1,times
!call psb_spmm(done,a,xv,dzero,bv,desc_a,info,'n')
call ap%spmm(done,xv,dzero,bv,info,trans='n')
end do
call psb_barrier(ictxt)
t2 = psb_wtime() - t1
call psb_amx(ictxt,t2)
! FIXME: cache flush needed here
call psb_barrier(ictxt)
tt1 = psb_wtime()
do i=1,times
call psb_spmm(done,a,xv,dzero,vtst,desc_a,info,trans='n')
end do
call psb_barrier(ictxt)
tt2 = psb_wtime() - tt1
call psb_amx(ictxt,tt2)
call psb_geaxpby(done,bv,-done,vtst,desc_a,info)
eps2 = psb_genrm2(vtst,desc_a,info)
epsi = psb_geamax(vtst,desc_a,info)
call psb_amx(ictxt,t2)
nr = desc_a%get_global_rows()
annz = a%get_nzeros()
amatsize = a%sizeof()
descsize = psb_sizeof(desc_a)
call psb_sum(ictxt,annz)
call psb_sum(ictxt,amatsize)
call psb_sum(ictxt,descsize)
if (iam == psb_root_) then
flops = 1.d1*2*annz
tflops=flops
write(psb_out_unit,'("Matrix: ell1 ",i0)') idim
write(psb_out_unit,'("Test on : ",i20," processors")') np
write(psb_out_unit,'("Size of matrix : ",i20," ")') nr
write(psb_out_unit,'("Number of nonzeros : ",i20," ")') annz
write(psb_out_unit,'("Memory occupation : ",i20," ")') amatsize
write(psb_out_unit,'("Number of flops (",i0," prod) : ",F20.0," ")') times,flops
flops = flops / (t2)
tflops = tflops / (tt2)
write(psb_out_unit,'("Time for ",i0," products (s) : ",F20.3)')times, t2
write(psb_out_unit,'("Time per product (ms) : ",F20.3)') t2*1.d3/(1.d0*times)
write(psb_out_unit,'("MFLOPS : ",F20.3)') flops/1.d6
write(psb_out_unit,'("Time for ",i0," products (s) (trans.): ",F20.3)') times,tt2
write(psb_out_unit,'("Time per product (ms) (trans.): ",F20.3)') tt2*1.d3/(1.d0*times)
write(psb_out_unit,'("MFLOPS (trans.): ",F20.3)') tflops/1.d6
write(psb_out_unit,'("DIFF nrm2 (trans.): ",F20.3)') eps2
write(psb_out_unit,'("DIFF nrmi (trans.): ",F20.3)') epsi
!
! This computation is valid for CSR
!
nbytes = nr*(2*psb_sizeof_dp + psb_sizeof_int)+&
& annz*(psb_sizeof_dp + psb_sizeof_int)
bdwdth = times*nbytes/(t2*1.d6)
write(psb_out_unit,*)
write(psb_out_unit,'("MBYTES/S : ",F20.3)') bdwdth
bdwdth = times*nbytes/(tt2*1.d6)
write(psb_out_unit,'("MBYTES/S (trans): ",F20.3)') bdwdth
write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt()
write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize
end if
!
! cleanup storage and exit
!
call psb_gefree(bv,desc_a,info)
call psb_gefree(xv,desc_a,info)
call psb_spfree(a,desc_a,info)
call psb_cdfree(desc_a,info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='free routine'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
9999 continue
if(info /= psb_success_) then
call psb_error(ictxt)
end if
call psb_exit(ictxt)
stop
contains
!
! get iteration parameters from standard input
!
subroutine get_parms(ictxt,afmt,idim)
integer(psb_ipk_) :: ictxt
character(len=*) :: afmt
integer(psb_ipk_) :: idim
integer(psb_ipk_) :: np, iam
integer(psb_ipk_) :: intbuf(10), ip
call psb_info(ictxt, iam, np)
if (iam == 0) then
read(psb_inp_unit,*) afmt
read(psb_inp_unit,*) idim
endif
call psb_bcast(ictxt,afmt)
call psb_bcast(ictxt,idim)
if (iam == 0) then
write(psb_out_unit,'("Testing matrix : ell1")')
write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim
write(psb_out_unit,'("Number of processors : ",i0)')np
write(psb_out_unit,'("Data distribution : BLOCK")')
write(psb_out_unit,'(" ")')
end if
return
end subroutine get_parms
!
! print an error message
!
subroutine pr_usage(iout)
integer(psb_ipk_) :: iout
write(iout,*)'incorrect parameter(s) found'
write(iout,*)' usage: pde90 methd prec dim &
&[istop itmax itrace]'
write(iout,*)' where:'
write(iout,*)' methd: cgstab cgs rgmres bicgstabl'
write(iout,*)' prec : bjac diag none'
write(iout,*)' dim number of points along each axis'
write(iout,*)' the size of the resulting linear '
write(iout,*)' system is dim**3'
write(iout,*)' istop stopping criterion 1, 2 '
write(iout,*)' itmax maximum number of iterations [500] '
write(iout,*)' itrace <=0 (no tracing, default) or '
write(iout,*)' >= 1 do tracing every itrace'
write(iout,*)' iterations '
end subroutine pr_usage
!
! functions parametrizing the differential equation
!
function b1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b1
real(psb_dpk_), intent(in) :: x,y,z
b1=1.d0/sqrt(3.d0)
end function b1
function b2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b2
real(psb_dpk_), intent(in) :: x,y,z
b2=1.d0/sqrt(3.d0)
end function b2
function b3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: b3
real(psb_dpk_), intent(in) :: x,y,z
b3=1.d0/sqrt(3.d0)
end function b3
function c(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: c
real(psb_dpk_), intent(in) :: x,y,z
c=0.d0
end function c
function a1(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a1
real(psb_dpk_), intent(in) :: x,y,z
a1=1.d0/80
end function a1
function a2(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a2
real(psb_dpk_), intent(in) :: x,y,z
a2=1.d0/80
end function a2
function a3(x,y,z)
use psb_base_mod, only : psb_dpk_
real(psb_dpk_) :: a3
real(psb_dpk_), intent(in) :: x,y,z
a3=1.d0/80
end function a3
function g(x,y,z)
use psb_base_mod, only : psb_dpk_, done
real(psb_dpk_) :: g
real(psb_dpk_), intent(in) :: x,y,z
g = dzero
if (x == done) then
g = done
else if (x == dzero) then
g = exp(y**2-z**2)
end if
end function g
end program pgenpmv

@ -2,7 +2,7 @@
BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES
BJAC Preconditioner NONE DIAG BJAC
CSR Storage format for matrix A: CSR COO JAD
100 Domain size (acutal system is this**3)
040 Domain size (acutal system is this**3)
2 Stopping criterion
1000 MAXIT
01 ITRACE

@ -55,7 +55,7 @@ subroutine psb_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
implicit none
procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
class(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
type(psb_desc_type) :: desc_a
integer(psb_ipk_) :: ictxt, info
@ -370,7 +370,7 @@ subroutine psb_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
implicit none
procedure(d_func_2d) :: b1,b2,c,a1,a2,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
class(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
type(psb_desc_type) :: desc_a
integer(psb_ipk_) :: ictxt, info
@ -391,8 +391,8 @@ subroutine psb_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
integer(psb_ipk_) :: np, iam, nr, nt
integer(psb_ipk_) :: icoeff
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_), allocatable :: irow(:),icol(:),myidx(:)
real(psb_dpk_), allocatable :: val(:)
! deltah dimension of each grid cell
! deltat discretization time
real(psb_dpk_) :: deltah, sqdeltah, deltah2

@ -66,7 +66,7 @@ module psb_d_genpde_mod
implicit none
procedure(d_func_3d) :: a1,a2,a3,c,b1,b2,b3,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
class(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
type(psb_desc_type) :: desc_a
integer(psb_ipk_) :: ictxt, info
@ -110,7 +110,7 @@ module psb_d_genpde_mod
implicit none
procedure(d_func_2d) :: a1,a2,c,b1,b2,g
integer(psb_ipk_) :: idim
type(psb_dspmat_type) :: a
class(psb_dspmat_type) :: a
type(psb_d_vect_type) :: xv,bv
type(psb_desc_type) :: desc_a
integer(psb_ipk_) :: ictxt, info

Loading…
Cancel
Save