From 39ed3311e750b15a4082fa840da736f7150ec40d Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 8 Oct 2013 15:26:55 +0000 Subject: [PATCH] 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. --- base/modules/Makefile | 3 + base/modules/psb_d_parmat_mod.f90 | 1095 +++++++++++++++++++++++++++++ base/modules/psb_d_psblas_mod.f90 | 4 +- base/modules/psb_d_tools_mod.f90 | 6 +- base/psblas/psb_dspmm.f90 | 206 +++++- base/tools/psb_dspalloc.f90 | 2 +- base/tools/psb_dspasb.f90 | 2 +- base/tools/psb_dspfree.f90 | 2 +- krylov/psb_dcgstab.f90 | 2 +- krylov/psb_dkrylov.f90 | 4 +- krylov/psb_krylov_mod.f90 | 2 +- test/kernel/Makefile | 6 +- test/kernel/pgenpmv.f90 | 321 +++++++++ test/kernel/runs/spmv.inp | 2 +- test/pargen/runs/ppde.inp | 2 +- util/psb_d_genpde_impl.f90 | 8 +- util/psb_d_genpde_mod.f90 | 4 +- 17 files changed, 1647 insertions(+), 24 deletions(-) create mode 100644 base/modules/psb_d_parmat_mod.f90 create mode 100644 test/kernel/pgenpmv.f90 diff --git a/base/modules/Makefile b/base/modules/Makefile index f4871876..454d8062 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -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) diff --git a/base/modules/psb_d_parmat_mod.f90 b/base/modules/psb_d_parmat_mod.f90 new file mode 100644 index 00000000..eccde436 --- /dev/null +++ b/base/modules/psb_d_parmat_mod.f90 @@ -0,0 +1,1095 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012, 2013 +!!$ 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. +!!$ +!!$ +! +! package: psb_d_mat_mod +! +! This module contains the definition of the psb_d_sparse type which +! is a generic container for a sparse matrix and it is mostly meant to +! provide a mean of switching, at run-time, among different formats, +! potentially unknown at the library compile-time by adding a layer of +! indirection. This type encapsulates the psb_d_base_sparse_mat class +! inside another class which is the one visible to the user. +! Most methods of the psb_d_mat_mod simply call the methods of the +! encapsulated class. +! The exceptions are mainly cscnv and cp_from/cp_to; these provide +! the functionalities to have the encapsulated class change its +! type dynamically, and to extract/input an inner object. +! +! A sparse matric has a state corresponding to its progression +! through the application life. +! In particular, computational methods can only be invoked when +! the matrix is in the ASSEMBLED state, whereas the other states are +! dedicated to operations on the internal matrix data. +! A sparse matrix can move between states according to the +! following state transition table. Associated with these states are +! the possible dynamic types of the inner matrix object. +! Only COO matrices can ever be in the BUILD state, whereas +! the ASSEMBLED and UPDATE state can be entered by any class. +! +! In Out Method +!| ---------------------------------- +!| Null Build csall +!| Build Build csput +!| Build Assembled cscnv +!| Assembled Assembled cscnv +!| Assembled Update reinit +!| Update Update csput +!| Update Assembled cscnv +!| * unchanged reall +!| Assembled Null free +! + + +module psb_d_parmat_mod + + use psb_d_mat_mod + use psb_desc_mod + + type :: psb_d_parmat_type + + type(psb_dspmat_type), allocatable :: a_d, a_nd + type(psb_desc_type), pointer :: desc + + contains + procedure, pass(a) :: create => psb_d_p_create + +!!$ ! Getters +!!$ procedure, pass(a) :: get_nrows => psb_d_p_get_nrows +!!$ procedure, pass(a) :: get_ncols => psb_d_p_get_ncols +!!$ procedure, pass(a) :: get_nzeros => psb_d_p_get_nzeros +!!$ procedure, pass(a) :: get_nz_row => psb_d_p_get_nz_row +!!$ procedure, pass(a) :: get_size => psb_d_p_get_size +!!$ procedure, pass(a) :: get_dupl => psb_d_p_get_dupl +!!$ procedure, pass(a) :: is_null => psb_d_p_is_null +!!$ procedure, pass(a) :: is_bld => psb_d_p_is_bld +!!$ procedure, pass(a) :: is_upd => psb_d_p_is_upd +!!$ procedure, pass(a) :: is_asb => psb_d_p_is_asb +!!$ procedure, pass(a) :: is_sorted => psb_d_p_is_sorted +!!$ procedure, pass(a) :: is_upper => psb_d_p_is_upper +!!$ procedure, pass(a) :: is_lower => psb_d_p_is_lower +!!$ procedure, pass(a) :: is_triangle => psb_d_p_is_triangle +!!$ procedure, pass(a) :: is_unit => psb_d_p_is_unit +!!$ procedure, pass(a) :: get_fmt => psb_d_p_get_fmt + procedure, pass(a) :: sizeof => psb_d_p_sizeof +!!$ +!!$ ! Setters +!!$ procedure, pass(a) :: set_nrows => psb_d_p_set_nrows +!!$ procedure, pass(a) :: set_ncols => psb_d_p_set_ncols +!!$ procedure, pass(a) :: set_dupl => psb_d_p_set_dupl +!!$ procedure, pass(a) :: set_null => psb_d_p_set_null +!!$ procedure, pass(a) :: set_bld => psb_d_p_set_bld +!!$ procedure, pass(a) :: set_upd => psb_d_p_set_upd +!!$ procedure, pass(a) :: set_asb => psb_d_p_set_asb +!!$ procedure, pass(a) :: set_sorted => psb_d_p_set_sorted +!!$ procedure, pass(a) :: set_upper => psb_d_p_set_upper +!!$ procedure, pass(a) :: set_lower => psb_d_p_set_lower +!!$ procedure, pass(a) :: set_triangle => psb_d_p_set_triangle +!!$ procedure, pass(a) :: set_unit => psb_d_p_set_unit +!!$ +!!$ ! Memory/data management +!!$ procedure, pass(a) :: csall => psb_d_p_csall +!!$ procedure, pass(a) :: free => psb_d_p_free +!!$ procedure, pass(a) :: trim => psb_d_p_trim +!!$ procedure, pass(a) :: csput => psb_d_p_csput +!!$ procedure, pass(a) :: csgetptn => psb_d_p_csgetptn +!!$ procedure, pass(a) :: csgetrow => psb_d_p_csgetrow +!!$ procedure, pass(a) :: csgetblk => psb_d_p_csgetblk +!!$ generic, public :: csget => csgetptn, csgetrow, csgetblk +!!$ procedure, pass(a) :: m_csclip => psb_d_p_csclip +!!$ procedure, pass(a) :: b_csclip => psb_d_p_b_csclip +!!$ generic, public :: csclip => b_csclip, m_csclip +!!$ procedure, pass(a) :: reall => psb_d_p_reallocate_nz +!!$ procedure, pass(a) :: get_neigh => psb_d_p_get_neigh +!!$ procedure, pass(a) :: reinit => psb_d_p_reinit +!!$ procedure, pass(a) :: print_i => psb_d_p_sparse_print +!!$ procedure, pass(a) :: print_n => psb_d_p_n_sparse_print +!!$ generic, public :: print => print_i, print_n +!!$ procedure, pass(a) :: mold => psb_d_p_mold +!!$ procedure, pass(a) :: transp_1mat => psb_d_p_transp_1mat +!!$ procedure, pass(a) :: transp_2mat => psb_d_p_transp_2mat +!!$ generic, public :: transp => transp_1mat, transp_2mat +!!$ procedure, pass(a) :: transc_1mat => psb_d_p_transc_1mat +!!$ procedure, pass(a) :: transc_2mat => psb_d_p_transc_2mat +!!$ generic, public :: transc => transc_1mat, transc_2mat +!!$ ! These are specific to this level of encapsulation. +!!$ procedure, pass(a) :: mv_from_b => psb_d_p_mv_from +!!$ generic, public :: mv_from => mv_from_b +!!$ procedure, pass(a) :: mv_to_b => psb_d_p_mv_to +!!$ generic, public :: mv_to => mv_to_b +!!$ procedure, pass(a) :: cp_from_b => psb_d_p_cp_from +!!$ generic, public :: cp_from => cp_from_b +!!$ procedure, pass(a) :: cp_to_b => psb_d_p_cp_to +!!$ generic, public :: cp_to => cp_to_b +!!$ procedure, pass(a) :: clip_d_ip => psb_d_p_clip_d_ip +!!$ procedure, pass(a) :: clip_d => psb_d_p_clip_d +!!$ generic, public :: clip_diag => clip_d_ip, clip_d +!!$ procedure, pass(a) :: cscnv_np => psb_d_p_cscnv +!!$ procedure, pass(a) :: cscnv_ip => psb_d_p_cscnv_ip +!!$ procedure, pass(a) :: cscnv_base => psb_d_p_cscnv_base +!!$ generic, public :: cscnv => cscnv_np, cscnv_ip, cscnv_base +!!$ procedure, pass(a) :: clone => psb_dspmat_clone +!!$ +!!$ ! Computational routines +!!$ procedure, pass(a) :: get_diag => psb_d_p_get_diag +!!$ procedure, pass(a) :: maxval => psb_d_p_maxval +!!$ procedure, pass(a) :: spnmi => psb_d_p_csnmi +!!$ procedure, pass(a) :: spnm1 => psb_d_p_csnm1 +!!$ procedure, pass(a) :: rowsum => psb_d_p_rowsum +!!$ procedure, pass(a) :: arwsum => psb_d_p_arwsum +!!$ procedure, pass(a) :: colsum => psb_d_p_colsum +!!$ procedure, pass(a) :: aclsum => psb_d_p_aclsum + procedure, pass(a) :: csmv_v => psb_d_p_csmv_vect +!!$ procedure, pass(a) :: csmv => psb_d_p_csmv +!!$ procedure, pass(a) :: csmm => psb_d_p_csmm + generic, public :: spmm => csmv_v +!!$ procedure, pass(a) :: scals => psb_d_p_scals +!!$ procedure, pass(a) :: scalv => psb_d_p_scal +!!$ generic, public :: scal => scals, scalv +!!$ procedure, pass(a) :: cssv_v => psb_d_p_cssv_vect +!!$ procedure, pass(a) :: cssv => psb_d_p_cssv +!!$ procedure, pass(a) :: cssm => psb_d_p_cssm +!!$ generic, public :: spsm => cssm, cssv, cssv_v + + end type psb_d_parmat_type + +!!$ private :: psb_d_get_nrows, psb_d_get_ncols, psb_d_get_nzeros, psb_d_get_size, & +!!$ & psb_d_get_dupl, psb_d_is_null, psb_d_is_bld, & +!!$ & psb_d_is_upd, psb_d_is_asb, psb_d_is_sorted, psb_d_is_upper, & +!!$ & psb_d_is_lower, psb_d_is_triangle, psb_d_get_nz_row + + interface psb_sizeof + module procedure psb_d_p_sizeof + end interface + + + ! == =================================== + ! + ! + ! + ! Setters + ! + ! + ! + ! + ! + ! + ! == =================================== + + +!!$ interface +!!$ subroutine psb_d_p_set_nrows(m,a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ integer(psb_ipk_), intent(in) :: m +!!$ end subroutine psb_d_p_set_nrows +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_ncols(n,a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ integer(psb_ipk_), intent(in) :: n +!!$ end subroutine psb_d_p_set_ncols +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_dupl(n,a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ integer(psb_ipk_), intent(in) :: n +!!$ end subroutine psb_d_p_set_dupl +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_null(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_set_null +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_bld(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_set_bld +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_upd(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_set_upd +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_asb(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_set_asb +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_sorted(a,val) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: val +!!$ end subroutine psb_d_p_set_sorted +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_triangle(a,val) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: val +!!$ end subroutine psb_d_p_set_triangle +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_unit(a,val) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: val +!!$ end subroutine psb_d_p_set_unit +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_lower(a,val) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: val +!!$ end subroutine psb_d_p_set_lower +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_set_upper(a,val) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: val +!!$ end subroutine psb_d_p_set_upper +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_sparse_print(iout,a,iv,head,ivr,ivc) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ integer(psb_ipk_), intent(in) :: iout +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_), intent(in), optional :: iv(:) +!!$ character(len=*), optional :: head +!!$ integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:) +!!$ end subroutine psb_d_p_sparse_print +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_n_sparse_print(fname,a,iv,head,ivr,ivc) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ character(len=*), intent(in) :: fname +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_), intent(in), optional :: iv(:) +!!$ character(len=*), optional :: head +!!$ integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:) +!!$ end subroutine psb_d_p_n_sparse_print +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_get_neigh(a,idx,neigh,n,info,lev) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_), intent(in) :: idx +!!$ integer(psb_ipk_), intent(out) :: n +!!$ integer(psb_ipk_), allocatable, intent(out) :: neigh(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), optional, intent(in) :: lev +!!$ end subroutine psb_d_p_get_neigh +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csall(nr,nc,a,info,nz) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(out) :: a +!!$ integer(psb_ipk_), intent(in) :: nr,nc +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: nz +!!$ end subroutine psb_d_p_csall +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_reallocate_nz(nz,a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ integer(psb_ipk_), intent(in) :: nz +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_reallocate_nz +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_free(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_free +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_trim(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_trim +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer(psb_ipk_), intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: gtl(:) +!!$ end subroutine psb_d_p_csput +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csgetptn(imin,imax,a,nz,ia,ja,info,& +!!$ & jmin,jmax,iren,append,nzin,rscale,cscale) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_), intent(in) :: imin,imax +!!$ integer(psb_ipk_), intent(out) :: nz +!!$ integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) +!!$ integer(psb_ipk_),intent(out) :: info +!!$ logical, intent(in), optional :: append +!!$ integer(psb_ipk_), intent(in), optional :: iren(:) +!!$ integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin +!!$ logical, intent(in), optional :: rscale,cscale +!!$ end subroutine psb_d_p_csgetptn +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csgetrow(imin,imax,a,nz,ia,ja,val,info,& +!!$ & jmin,jmax,iren,append,nzin,rscale,cscale) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_), intent(in) :: imin,imax +!!$ integer(psb_ipk_), intent(out) :: nz +!!$ integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) +!!$ real(psb_dpk_), allocatable, intent(inout) :: val(:) +!!$ integer(psb_ipk_),intent(out) :: info +!!$ logical, intent(in), optional :: append +!!$ integer(psb_ipk_), intent(in), optional :: iren(:) +!!$ integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin +!!$ logical, intent(in), optional :: rscale,cscale +!!$ end subroutine psb_d_p_csgetrow +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csgetblk(imin,imax,a,b,info,& +!!$ & jmin,jmax,iren,append,rscale,cscale) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dspmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_dspmat_type), intent(inout) :: b +!!$ integer(psb_ipk_), intent(in) :: imin,imax +!!$ integer(psb_ipk_),intent(out) :: info +!!$ logical, intent(in), optional :: append +!!$ integer(psb_ipk_), intent(in), optional :: iren(:) +!!$ integer(psb_ipk_), intent(in), optional :: jmin,jmax +!!$ logical, intent(in), optional :: rscale,cscale +!!$ end subroutine psb_d_p_csgetblk +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_csclip(a,b,info,& +!!$ & imin,imax,jmin,jmax,rscale,cscale) +!!$ import :: psb_ipk_, psb_dspmat_type, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_dspmat_type), intent(out) :: b +!!$ integer(psb_ipk_),intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax +!!$ logical, intent(in), optional :: rscale,cscale +!!$ end subroutine psb_d_p_csclip +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_b_csclip(a,b,info,& +!!$ & imin,imax,jmin,jmax,rscale,cscale) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_coo_sparse_mat +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ type(psb_d_coo_sparse_mat), intent(out) :: b +!!$ integer(psb_ipk_),intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax +!!$ logical, intent(in), optional :: rscale,cscale +!!$ end subroutine psb_d_p_b_csclip +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_mold(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ class(psb_d_base_sparse_mat), allocatable, intent(out) :: b +!!$ end subroutine psb_d_p_mold +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_transp_1mat(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_transp_1mat +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_transp_2mat(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_parmat_type), intent(out) :: b +!!$ end subroutine psb_d_p_transp_2mat +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_transc_1mat(a) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ end subroutine psb_d_p_transc_1mat +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_transc_2mat(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_parmat_type), intent(out) :: b +!!$ end subroutine psb_d_p_transc_2mat +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_reinit(a,clear) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ logical, intent(in), optional :: clear +!!$ end subroutine psb_d_p_reinit +!!$ +!!$ end interface +!!$ +!!$ +!!$ ! +!!$ ! These methods are specific to the outer SPMAT_TYPE level, since +!!$ ! they tamper with the inner BASE_SPARSE_MAT object. +!!$ ! +!!$ ! +!!$ +!!$ ! +!!$ ! CSCNV: switches to a different internal derived type. +!!$ ! 3 versions: copying to target +!!$ ! copying to a base_sparse_mat object. +!!$ ! in place +!!$ ! +!!$ ! +!!$ interface +!!$ subroutine psb_d_p_cscnv(a,b,info,type,mold,upd,dupl) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_parmat_type), intent(out) :: b +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_),optional, intent(in) :: dupl, upd +!!$ character(len=*), optional, intent(in) :: type +!!$ class(psb_d_base_sparse_mat), intent(in), optional :: mold +!!$ end subroutine psb_d_p_cscnv +!!$ end interface +!!$ +!!$ +!!$ interface +!!$ subroutine psb_d_p_cscnv_ip(a,iinfo,type,mold,dupl) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ integer(psb_ipk_), intent(out) :: iinfo +!!$ integer(psb_ipk_),optional, intent(in) :: dupl +!!$ character(len=*), optional, intent(in) :: type +!!$ class(psb_d_base_sparse_mat), intent(in), optional :: mold +!!$ end subroutine psb_d_p_cscnv_ip +!!$ end interface +!!$ +!!$ +!!$ interface +!!$ subroutine psb_d_p_cscnv_base(a,b,info,dupl) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_base_sparse_mat), intent(out) :: b +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_),optional, intent(in) :: dupl +!!$ end subroutine psb_d_p_cscnv_base +!!$ end interface +!!$ +!!$ ! +!!$ ! Produce a version of the matrix with diagonal cut +!!$ ! out; passes through a COO buffer. +!!$ ! +!!$ interface +!!$ subroutine psb_d_p_clip_d(a,b,info) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_parmat_type), intent(out) :: b +!!$ integer(psb_ipk_),intent(out) :: info +!!$ end subroutine psb_d_p_clip_d +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_clip_d_ip(a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ integer(psb_ipk_),intent(out) :: info +!!$ end subroutine psb_d_p_clip_d_ip +!!$ end interface +!!$ +!!$ ! +!!$ ! These four interfaces cut through the +!!$ ! encapsulation between spmat_type and base_sparse_mat. +!!$ ! +!!$ interface +!!$ subroutine psb_d_p_mv_from(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(out) :: a +!!$ class(psb_d_base_sparse_mat), intent(inout) :: b +!!$ end subroutine psb_d_p_mv_from +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_cp_from(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(out) :: a +!!$ class(psb_d_base_sparse_mat), intent(in) :: b +!!$ end subroutine psb_d_p_cp_from +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_mv_to(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ class(psb_d_base_sparse_mat), intent(out) :: b +!!$ end subroutine psb_d_p_mv_to +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_cp_to(a,b) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_, psb_d_base_sparse_mat +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ class(psb_d_base_sparse_mat), intent(out) :: b +!!$ end subroutine psb_d_p_cp_to +!!$ end interface +!!$ +!!$ ! +!!$ ! Transfer the internal allocation to the target. +!!$ ! +!!$ interface psb_move_alloc +!!$ subroutine psb_d_p_parmat_type_move(a,b,info) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ class(psb_d_parmat_type), intent(out) :: b +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_parmat_type_move +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_pspmat_clone(a,b,info) +!!$ import :: psb_ipk_, psb_d_parmat_type +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ class(psb_d_parmat_type), intent(inout) :: b +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_pspmat_clone +!!$ end interface +!!$ + + + + ! == =================================== + ! + ! + ! + ! Computational routines + ! + ! + ! + ! + ! + ! + ! == =================================== + + interface + subroutine psb_d_p_csmv_vect(alpha,a,x,beta,y,info,trans) + use psb_d_vect_mod, only : psb_d_vect_type + import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ + 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 + end subroutine psb_d_p_csmv_vect + end interface + +!!$ interface psb_csmm +!!$ subroutine psb_d_p_csmm(alpha,a,x,beta,y,info,trans) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) +!!$ real(psb_dpk_), intent(inout) :: y(:,:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character, optional, intent(in) :: trans +!!$ end subroutine psb_d_p_csmm +!!$ subroutine psb_d_p_csmv(alpha,a,x,beta,y,info,trans) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(in) :: alpha, beta, x(:) +!!$ real(psb_dpk_), intent(inout) :: y(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character, optional, intent(in) :: trans +!!$ end subroutine psb_d_p_csmv +!!$ subroutine psb_d_p_csmv_vect(alpha,a,x,beta,y,info,trans) +!!$ use psb_d_vect_mod, only : psb_d_vect_type +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ 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 +!!$ end subroutine psb_d_p_csmv_vect +!!$ end interface +!!$ +!!$ interface psb_cssm +!!$ subroutine psb_d_p_cssm(alpha,a,x,beta,y,info,trans,scale,d) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) +!!$ real(psb_dpk_), intent(inout) :: y(:,:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character, optional, intent(in) :: trans, scale +!!$ real(psb_dpk_), intent(in), optional :: d(:) +!!$ end subroutine psb_d_p_cssm +!!$ subroutine psb_d_p_cssv(alpha,a,x,beta,y,info,trans,scale,d) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(in) :: alpha, beta, x(:) +!!$ real(psb_dpk_), intent(inout) :: y(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character, optional, intent(in) :: trans, scale +!!$ real(psb_dpk_), intent(in), optional :: d(:) +!!$ end subroutine psb_d_p_cssv +!!$ subroutine psb_d_p_cssv_vect(alpha,a,x,beta,y,info,trans,scale,d) +!!$ use psb_d_vect_mod, only : psb_d_vect_type +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ 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, scale +!!$ type(psb_d_vect_type), optional, intent(inout) :: d +!!$ end subroutine psb_d_p_cssv_vect +!!$ end interface +!!$ +!!$ interface +!!$ function psb_d_maxval(a) result(res) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_) :: res +!!$ end function psb_d_maxval +!!$ end interface +!!$ +!!$ interface +!!$ function psb_d_csnmi(a) result(res) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_) :: res +!!$ end function psb_d_csnmi +!!$ end interface +!!$ +!!$ interface +!!$ function psb_d_csnm1(a) result(res) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_) :: res +!!$ end function psb_d_csnm1 +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_rowsum(d,a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(out) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_rowsum +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_arwsum(d,a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(out) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_arwsum +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_colsum(d,a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(out) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_colsum +!!$ end interface +!!$ +!!$ interface +!!$ subroutine psb_d_p_aclsum(d,a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(out) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_aclsum +!!$ end interface +!!$ +!!$ +!!$ interface +!!$ subroutine psb_d_p_get_diag(a,d,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ real(psb_dpk_), intent(out) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_get_diag +!!$ end interface +!!$ +!!$ interface psb_scal +!!$ subroutine psb_d_p_scal(d,a,info,side) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ real(psb_dpk_), intent(in) :: d(:) +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character, intent(in), optional :: side +!!$ end subroutine psb_d_p_scal +!!$ subroutine psb_d_p_scals(d,a,info) +!!$ import :: psb_ipk_, psb_d_parmat_type, psb_dpk_ +!!$ class(psb_d_parmat_type), intent(inout) :: a +!!$ real(psb_dpk_), intent(in) :: d +!!$ integer(psb_ipk_), intent(out) :: info +!!$ end subroutine psb_d_p_scals +!!$ end interface + + +contains + + + ! == =================================== + ! + ! + ! + ! Getters + ! + ! + ! + ! + ! + ! == =================================== + + + function psb_d_p_sizeof(a) result(res) + implicit none + class(psb_d_parmat_type), intent(in) :: a + integer(psb_long_int_k_) :: res + + res = 0 + if (allocated(a%a_d)) & + & res = res + a%a_d%sizeof() + if (allocated(a%a_nd)) & + & res = res + a%a_nd%sizeof() + + end function psb_d_p_sizeof + + +!!$ function psb_d_p_get_fmt(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ character(len=5) :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_fmt() +!!$ else +!!$ res = 'NULL' +!!$ end if +!!$ +!!$ end function psb_d_p_get_fmt +!!$ +!!$ +!!$ function psb_d_p_get_dupl(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_dupl() +!!$ else +!!$ res = psb_invalid_ +!!$ end if +!!$ end function psb_d_p_get_dupl +!!$ +!!$ function psb_d_p_get_nrows(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_nrows() +!!$ else +!!$ res = 0 +!!$ end if +!!$ +!!$ end function psb_d_p_get_nrows +!!$ +!!$ function psb_d_p_get_ncols(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_ncols() +!!$ else +!!$ res = 0 +!!$ end if +!!$ +!!$ end function psb_d_p_get_ncols +!!$ +!!$ function psb_d_p_is_triangle(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_triangle() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_triangle +!!$ +!!$ function psb_d_p_is_unit(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_unit() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_unit +!!$ +!!$ function psb_d_p_is_upper(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_upper() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_upper +!!$ +!!$ function psb_d_p_is_lower(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = .not. a%a%is_upper() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_lower +!!$ +!!$ function psb_d_p_is_null(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_null() +!!$ else +!!$ res = .true. +!!$ end if +!!$ +!!$ end function psb_d_p_is_null +!!$ +!!$ function psb_d_p_is_bld(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_bld() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_bld +!!$ +!!$ function psb_d_p_is_upd(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_upd() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_upd +!!$ +!!$ function psb_d_p_is_asb(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_asb() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_asb +!!$ +!!$ function psb_d_p_is_sorted(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ logical :: res +!!$ +!!$ if (allocated(a%a)) then +!!$ res = a%a%is_sorted() +!!$ else +!!$ res = .false. +!!$ end if +!!$ +!!$ end function psb_d_p_is_sorted +!!$ +!!$ +!!$ +!!$ function psb_d_p_get_nzeros(a) result(res) +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = 0 +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_nzeros() +!!$ end if +!!$ +!!$ end function psb_d_p_get_nzeros +!!$ +!!$ function psb_d_p_get_size(a) result(res) +!!$ +!!$ implicit none +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ +!!$ res = 0 +!!$ if (allocated(a%a)) then +!!$ res = a%a%get_size() +!!$ end if +!!$ +!!$ end function psb_d_p_get_size +!!$ +!!$ +!!$ function psb_d_p_get_nz_row(idx,a) result(res) +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: idx +!!$ class(psb_d_parmat_type), intent(in) :: a +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = 0 +!!$ +!!$ if (allocated(a%a)) res = a%a%get_nz_row(idx) +!!$ +!!$ end function psb_d_p_get_nz_row + + subroutine psb_d_p_create(a,as,desc,info) + use psb_error_mod + implicit none + class(psb_d_parmat_type), intent(out) :: a + type(psb_dspmat_type), intent(in) :: as + type(psb_desc_type), intent(in), target :: desc + integer(psb_ipk_), intent(out) :: info + character(len=20), parameter :: name='psb_d_p_create' + integer(psb_ipk_) :: err_act, nr, nc + integer(psb_ipk_) :: ierr(5) + + call psb_erractionsave(err_act) + info = psb_success_ + + a%desc => desc + allocate(a%a_d, a%a_nd,stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + nr = desc%get_local_rows() + nc = desc%get_local_cols() + if ((nr /= as%get_nrows()).or.(nc /= as%get_ncols())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + if (info == 0) call as%csclip(a%a_d,info,imax=nr,jmax=nr) + if (info == 0) call a%a_d%cscnv(info,mold=as%a) + if (info == 0) call as%csclip(a%a_nd,info,imax=nr,jmin=nr+1,jmax=nc) + if (info == 0) call a%a_nd%cscnv(info,mold=as%a) + if (info /= 0) then + info = psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 + endif + + 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 + + + end subroutine psb_d_p_create + + + +end module psb_d_parmat_mod diff --git a/base/modules/psb_d_psblas_mod.f90 b/base/modules/psb_d_psblas_mod.f90 index 70f0fdac..4e021545 100644 --- a/base/modules/psb_d_psblas_mod.f90 +++ b/base/modules/psb_d_psblas_mod.f90 @@ -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 diff --git a/base/modules/psb_d_tools_mod.f90 b/base/modules/psb_d_tools_mod.f90 index 3522bbb7..14bbd5a7 100644 --- a/base/modules/psb_d_tools_mod.f90 +++ b/base/modules/psb_d_tools_mod.f90 @@ -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 diff --git a/base/psblas/psb_dspmm.f90 b/base/psblas/psb_dspmm.f90 index ec34195e..582d488d 100644 --- a/base/psblas/psb_dspmm.f90 +++ b/base/psblas/psb_dspmm.f90 @@ -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 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 + diff --git a/base/tools/psb_dspalloc.f90 b/base/tools/psb_dspalloc.f90 index 6cf9e197..93b347b8 100644 --- a/base/tools/psb_dspalloc.f90 +++ b/base/tools/psb_dspalloc.f90 @@ -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 diff --git a/base/tools/psb_dspasb.f90 b/base/tools/psb_dspasb.f90 index ffa09f12..edee4d7a 100644 --- a/base/tools/psb_dspasb.f90 +++ b/base/tools/psb_dspasb.f90 @@ -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 diff --git a/base/tools/psb_dspfree.f90 b/base/tools/psb_dspfree.f90 index 1352a72f..be9bc931 100644 --- a/base/tools/psb_dspfree.f90 +++ b/base/tools/psb_dspfree.f90 @@ -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 diff --git a/krylov/psb_dcgstab.f90 b/krylov/psb_dcgstab.f90 index 03bc8c50..f22e38f8 100644 --- a/krylov/psb_dcgstab.f90 +++ b/krylov/psb_dcgstab.f90 @@ -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 diff --git a/krylov/psb_dkrylov.f90 b/krylov/psb_dkrylov.f90 index 5e75d9a3..23baaba8 100644 --- a/krylov/psb_dkrylov.f90 +++ b/krylov/psb_dkrylov.f90 @@ -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 diff --git a/krylov/psb_krylov_mod.f90 b/krylov/psb_krylov_mod.f90 index 78c474ba..1f709d4d 100644 --- a/krylov/psb_krylov_mod.f90 +++ b/krylov/psb_krylov_mod.f90 @@ -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 diff --git a/test/kernel/Makefile b/test/kernel/Makefile index 50b27000..0ecf332a 100644 --- a/test/kernel/Makefile +++ b/test/kernel/Makefile @@ -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) diff --git a/test/kernel/pgenpmv.f90 b/test/kernel/pgenpmv.f90 new file mode 100644 index 00000000..caf3dea6 --- /dev/null +++ b/test/kernel/pgenpmv.f90 @@ -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 diff --git a/test/kernel/runs/spmv.inp b/test/kernel/runs/spmv.inp index 3b4e9323..72025415 100644 --- a/test/kernel/runs/spmv.inp +++ b/test/kernel/runs/spmv.inp @@ -1,3 +1,3 @@ CSR -10 +80 diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index ad9aac03..e86827e4 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -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 diff --git a/util/psb_d_genpde_impl.f90 b/util/psb_d_genpde_impl.f90 index 989ebb40..4c3ea2a3 100644 --- a/util/psb_d_genpde_impl.f90 +++ b/util/psb_d_genpde_impl.f90 @@ -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 diff --git a/util/psb_d_genpde_mod.f90 b/util/psb_d_genpde_mod.f90 index 1afafec7..716f7e7f 100644 --- a/util/psb_d_genpde_mod.f90 +++ b/util/psb_d_genpde_mod.f90 @@ -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