From 40755d62791723521f08f187eeb9a39f853eab57 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sun, 8 Sep 2019 16:52:20 +0100 Subject: [PATCH] New method par_csr_spspmm. --- base/modules/tools/psb_c_tools_mod.f90 | 16 +- base/modules/tools/psb_d_tools_mod.f90 | 16 +- base/modules/tools/psb_s_tools_mod.f90 | 16 +- base/modules/tools/psb_z_tools_mod.f90 | 16 +- base/tools/Makefile | 3 +- base/tools/psb_c_par_csr_spspmm.f90 | 249 +++++++++++++++++++++++++ base/tools/psb_d_par_csr_spspmm.f90 | 249 +++++++++++++++++++++++++ base/tools/psb_s_par_csr_spspmm.f90 | 249 +++++++++++++++++++++++++ base/tools/psb_z_par_csr_spspmm.f90 | 249 +++++++++++++++++++++++++ configure.ac | 5 - 10 files changed, 1058 insertions(+), 10 deletions(-) create mode 100644 base/tools/psb_c_par_csr_spspmm.f90 create mode 100644 base/tools/psb_d_par_csr_spspmm.f90 create mode 100644 base/tools/psb_s_par_csr_spspmm.f90 create mode 100644 base/tools/psb_z_par_csr_spspmm.f90 diff --git a/base/modules/tools/psb_c_tools_mod.f90 b/base/modules/tools/psb_c_tools_mod.f90 index e31ef96a..c190521b 100644 --- a/base/modules/tools/psb_c_tools_mod.f90 +++ b/base/modules/tools/psb_c_tools_mod.f90 @@ -32,7 +32,7 @@ Module psb_c_tools_mod use psb_desc_mod, only : psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ use psb_c_vect_mod, only : psb_c_base_vect_type, psb_c_vect_type - use psb_c_mat_mod, only : psb_cspmat_type, psb_lcspmat_type, psb_c_base_sparse_mat + use psb_c_mat_mod, only : psb_cspmat_type, psb_lcspmat_type, psb_c_base_sparse_mat, psb_lc_csr_sparse_mat use psb_l_vect_mod, only : psb_l_vect_type use psb_c_multivect_mod, only : psb_c_base_multivect_type, psb_c_multivect_type @@ -319,5 +319,19 @@ Module psb_c_tools_mod logical, intent(in), optional :: clear end subroutine psb_csprn end interface + + interface psb_par_spspmm + subroutine psb_lc_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + import :: psb_lc_csr_sparse_mat, psb_desc_type, psb_ipk_ + Implicit None + type(psb_lc_csr_sparse_mat),intent(in) :: acsr + type(psb_lc_csr_sparse_mat),intent(inout) :: bcsr + type(psb_lc_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + End Subroutine psb_lc_par_csr_spspmm + end interface psb_par_spspmm end module psb_c_tools_mod diff --git a/base/modules/tools/psb_d_tools_mod.f90 b/base/modules/tools/psb_d_tools_mod.f90 index 971fc6b9..ce3b4210 100644 --- a/base/modules/tools/psb_d_tools_mod.f90 +++ b/base/modules/tools/psb_d_tools_mod.f90 @@ -32,7 +32,7 @@ Module psb_d_tools_mod use psb_desc_mod, only : psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ use psb_d_vect_mod, only : psb_d_base_vect_type, psb_d_vect_type - use psb_d_mat_mod, only : psb_dspmat_type, psb_ldspmat_type, psb_d_base_sparse_mat + use psb_d_mat_mod, only : psb_dspmat_type, psb_ldspmat_type, psb_d_base_sparse_mat, psb_ld_csr_sparse_mat use psb_l_vect_mod, only : psb_l_vect_type use psb_d_multivect_mod, only : psb_d_base_multivect_type, psb_d_multivect_type @@ -319,5 +319,19 @@ Module psb_d_tools_mod logical, intent(in), optional :: clear end subroutine psb_dsprn end interface + + interface psb_par_spspmm + subroutine psb_ld_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + import :: psb_ld_csr_sparse_mat, psb_desc_type, psb_ipk_ + Implicit None + type(psb_ld_csr_sparse_mat),intent(in) :: acsr + type(psb_ld_csr_sparse_mat),intent(inout) :: bcsr + type(psb_ld_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + End Subroutine psb_ld_par_csr_spspmm + end interface psb_par_spspmm end module psb_d_tools_mod diff --git a/base/modules/tools/psb_s_tools_mod.f90 b/base/modules/tools/psb_s_tools_mod.f90 index 155a39d0..7bc0e791 100644 --- a/base/modules/tools/psb_s_tools_mod.f90 +++ b/base/modules/tools/psb_s_tools_mod.f90 @@ -32,7 +32,7 @@ Module psb_s_tools_mod use psb_desc_mod, only : psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_ use psb_s_vect_mod, only : psb_s_base_vect_type, psb_s_vect_type - use psb_s_mat_mod, only : psb_sspmat_type, psb_lsspmat_type, psb_s_base_sparse_mat + use psb_s_mat_mod, only : psb_sspmat_type, psb_lsspmat_type, psb_s_base_sparse_mat, psb_ls_csr_sparse_mat use psb_l_vect_mod, only : psb_l_vect_type use psb_s_multivect_mod, only : psb_s_base_multivect_type, psb_s_multivect_type @@ -319,5 +319,19 @@ Module psb_s_tools_mod logical, intent(in), optional :: clear end subroutine psb_ssprn end interface + + interface psb_par_spspmm + subroutine psb_ls_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + import :: psb_ls_csr_sparse_mat, psb_desc_type, psb_ipk_ + Implicit None + type(psb_ls_csr_sparse_mat),intent(in) :: acsr + type(psb_ls_csr_sparse_mat),intent(inout) :: bcsr + type(psb_ls_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + End Subroutine psb_ls_par_csr_spspmm + end interface psb_par_spspmm end module psb_s_tools_mod diff --git a/base/modules/tools/psb_z_tools_mod.f90 b/base/modules/tools/psb_z_tools_mod.f90 index 565ff8d1..1b1ce21e 100644 --- a/base/modules/tools/psb_z_tools_mod.f90 +++ b/base/modules/tools/psb_z_tools_mod.f90 @@ -32,7 +32,7 @@ Module psb_z_tools_mod use psb_desc_mod, only : psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ use psb_z_vect_mod, only : psb_z_base_vect_type, psb_z_vect_type - use psb_z_mat_mod, only : psb_zspmat_type, psb_lzspmat_type, psb_z_base_sparse_mat + use psb_z_mat_mod, only : psb_zspmat_type, psb_lzspmat_type, psb_z_base_sparse_mat, psb_lz_csr_sparse_mat use psb_l_vect_mod, only : psb_l_vect_type use psb_z_multivect_mod, only : psb_z_base_multivect_type, psb_z_multivect_type @@ -319,5 +319,19 @@ Module psb_z_tools_mod logical, intent(in), optional :: clear end subroutine psb_zsprn end interface + + interface psb_par_spspmm + subroutine psb_lz_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + import :: psb_lz_csr_sparse_mat, psb_desc_type, psb_ipk_ + Implicit None + type(psb_lz_csr_sparse_mat),intent(in) :: acsr + type(psb_lz_csr_sparse_mat),intent(inout) :: bcsr + type(psb_lz_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + End Subroutine psb_lz_par_csr_spspmm + end interface psb_par_spspmm end module psb_z_tools_mod diff --git a/base/tools/Makefile b/base/tools/Makefile index 962a7faa..c29bee88 100644 --- a/base/tools/Makefile +++ b/base/tools/Makefile @@ -24,7 +24,8 @@ FOBJS = psb_cdall.o psb_cdals.o psb_cdalv.o psb_cd_inloc.o psb_cdins.o psb_cdprt psb_zspins.o psb_zsprn.o \ psb_cspalloc.o psb_cspasb.o psb_cspfree.o\ psb_cspins.o psb_csprn.o psb_cd_set_bld.o \ - psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o + psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o \ + psb_s_par_csr_spspmm.o psb_d_par_csr_spspmm.o psb_c_par_csr_spspmm.o psb_z_par_csr_spspmm.o # psb_lallc.o psb_lasb.o psb_lfree.o psb_lins.o \ MPFOBJS = psb_icdasb.o psb_ssphalo.o psb_dsphalo.o psb_csphalo.o psb_zsphalo.o \ diff --git a/base/tools/psb_c_par_csr_spspmm.f90 b/base/tools/psb_c_par_csr_spspmm.f90 new file mode 100644 index 00000000..d835e27d --- /dev/null +++ b/base/tools/psb_c_par_csr_spspmm.f90 @@ -0,0 +1,249 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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: psb_c_par_csr_spspmm.f90 +! +! Subroutine: psb_c_par_csr_spspmm +! Version: complex +! +! This routine computes a parallel product of two sparse matrices +! +! C = A * B +! +! where all the matrices are stored in CSR. On input and output the matrices +! are stored with column indices in local numbering, but inermediate quantities +! are in global numbering because gathering the halo of B to multiply it +! by A implies a potential enlargement of the support. +! Also, B may have a column index space different from its row index space, +! which is obviously the same as the column space of A. +! +! +! Arguments: +! acsr - type(psb_c_csr_sparse_mat), input. +! The sparse matrix structure A +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the column space of A +! bcsr - type(psb_c_csr_sparse_mat), input/output. +! The sparse matrix structure B, gets row-extended on output +! ccsr - type(psb_c_csr_sparse_mat), output +! The sparse matrix structure C +! desc_c - type(psb_desc_type), input/output. +! The communication descriptor of the column space of B +! +! info - integer, output. +! Error code. +! +!!$Subroutine psb_c_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) +!!$ use psb_base_mod, psb_protect_name => psb_c_par_csr_spspmm +!!$ Implicit None +!!$ +!!$ type(psb_c_csr_sparse_mat),intent(in) :: acsr +!!$ type(psb_c_csr_sparse_mat),intent(inout) :: bcsr +!!$ type(psb_c_csr_sparse_mat),intent(out) :: ccsr +!!$ type(psb_desc_type),intent(in) :: desc_a +!!$ type(psb_desc_type),intent(inout) :: desc_c +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: data +!!$ ! ...local scalars.... +!!$ integer(psb_ipk_) :: ictxt, np,me +!!$ integer(psb_ipk_) :: ncol, nnz +!!$ type(psb_c_csr_sparse_mat) :: tcsr1 +!!$ logical :: update_desc_c +!!$ integer(psb_ipk_) :: debug_level, debug_unit, err_act +!!$ character(len=20) :: name, ch_err +!!$ +!!$ if(psb_get_errstatus() /= 0) return +!!$ info=psb_success_ +!!$ name='psb_c_p_csr_spspmm' +!!$ call psb_erractionsave(err_act) +!!$ if (psb_errstatus_fatal()) then +!!$ info = psb_err_internal_error_ ; goto 9999 +!!$ end if +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ ictxt = desc_a%get_context() +!!$ +!!$ call psb_info(ictxt, me, np) +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),': Start' +!!$ +!!$ update_desc_c = desc_c%is_bld() +!!$ +!!$ ! +!!$ ! This is a bit tricky. +!!$ ! DESC_A is the descriptor of (the columns of) A, and therefore +!!$ ! of the rows of B; the columns of B, in the intended usage, span +!!$ ! a different space for which we have DESC_C. +!!$ ! We are gathering the halo rows of B to multiply by A; +!!$ ! now, the columns of B would ideally be kept in +!!$ ! global numbering, so that we can call this repeatedly to accumulate +!!$ ! the product of multiple operators, and convert to local numbering +!!$ ! at the last possible moment. However, this would imply calling +!!$ ! the serial SPSPMM with a matrix B with the GLOBAL number of columns +!!$ ! and this could be very expensive in memory. The solution is to keep B +!!$ ! in local numbering, so that only columns really appearing count, but to +!!$ ! expand the descriptor when gathering the halo, because by performing +!!$ ! the products we are extending the support of the operator; hence +!!$ ! this routine is intended to be called with a temporary descriptor +!!$ ! DESC_C which is in the BUILD state, to allow for such expansion +!!$ ! across multiple products. +!!$ ! The caller will at some later point finalize the descriptor DESC_C. +!!$ ! +!!$ +!!$ ncol = desc_a%get_local_cols() +!!$ call psb_sphalo(bcsr,desc_a,tcsr1,info,& +!!$ & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) +!!$ nnz = tcsr1%get_nzeros() +!!$ if (update_desc_c) then +!!$ call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) +!!$ else +!!$ call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) +!!$ end if +!!$ if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1) +!!$ if (info == psb_success_) call tcsr1%free() +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ call bcsr%set_ncols(desc_c%get_local_cols()) +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm 3' +!!$ if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() +!!$ call psb_spspmm(acsr,bcsr,ccsr,info) +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(ictxt,err_act) +!!$ +!!$ return +!!$ +!!$End Subroutine psb_c_par_csr_spspmm + +Subroutine psb_lc_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + use psb_base_mod, psb_protect_name => psb_lc_par_csr_spspmm + Implicit None + + type(psb_lc_csr_sparse_mat),intent(in) :: acsr + type(psb_lc_csr_sparse_mat),intent(inout) :: bcsr + type(psb_lc_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + ! ...local scalars.... + integer(psb_ipk_) :: ictxt, np,me + integer(psb_lpk_) :: nacol, nccol, nnz + type(psb_lc_csr_sparse_mat) :: tcsr1 + logical :: update_desc_c + integer(psb_ipk_) :: debug_level, debug_unit, err_act + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='psb_lc_p_csr_spspmm' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_a%get_context() + + call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + update_desc_c = desc_c%is_bld() + + ! + ! This is a bit tricky. + ! DESC_A is the descriptor of (the columns of) A, and therefore + ! of the rows of B; the columns of B, in the intended usage, span + ! a different space for which we have DESC_C. + ! We are gathering the halo rows of B to multiply by A; + ! now, the columns of B would ideally be kept in + ! global numbering, so that we can call this repeatedly to accumulate + ! the product of multiple operators, and convert to local numbering + ! at the last possible moment. However, this would imply calling + ! the serial SPSPMM with a matrix B with the GLOBAL number of columns + ! and this could be very expensive in memory. The solution is to keep B + ! in local numbering, so that only columns really appearing count, but to + ! expand the descriptor when gathering the halo, because by performing + ! the products we are extending the support of the operator; hence + ! this routine is intended to be called with a temporary descriptor + ! DESC_C which is in the BUILD state, to allow for such expansion + ! across multiple products. + ! The caller will at some later point finalize the descriptor DESC_C. + ! + + nacol = desc_a%get_local_cols() + call psb_sphalo(bcsr,desc_a,tcsr1,info,& + & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) + nnz = tcsr1%get_nzeros() + if (update_desc_c) then + call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) + else + call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) + end if + if (info == psb_success_) call psb_rwextd(nacol,bcsr,info,b=tcsr1) + if (info == psb_success_) call tcsr1%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + nccol = desc_c%get_local_cols() + call bcsr%set_ncols(nccol) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm 3' + if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() + call psb_spspmm(acsr,bcsr,ccsr,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +End Subroutine psb_lc_par_csr_spspmm diff --git a/base/tools/psb_d_par_csr_spspmm.f90 b/base/tools/psb_d_par_csr_spspmm.f90 new file mode 100644 index 00000000..7edfbac8 --- /dev/null +++ b/base/tools/psb_d_par_csr_spspmm.f90 @@ -0,0 +1,249 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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: psb_d_par_csr_spspmm.f90 +! +! Subroutine: psb_d_par_csr_spspmm +! Version: real +! +! This routine computes a parallel product of two sparse matrices +! +! C = A * B +! +! where all the matrices are stored in CSR. On input and output the matrices +! are stored with column indices in local numbering, but inermediate quantities +! are in global numbering because gathering the halo of B to multiply it +! by A implies a potential enlargement of the support. +! Also, B may have a column index space different from its row index space, +! which is obviously the same as the column space of A. +! +! +! Arguments: +! acsr - type(psb_d_csr_sparse_mat), input. +! The sparse matrix structure A +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the column space of A +! bcsr - type(psb_d_csr_sparse_mat), input/output. +! The sparse matrix structure B, gets row-extended on output +! ccsr - type(psb_d_csr_sparse_mat), output +! The sparse matrix structure C +! desc_c - type(psb_desc_type), input/output. +! The communication descriptor of the column space of B +! +! info - integer, output. +! Error code. +! +!!$Subroutine psb_d_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) +!!$ use psb_base_mod, psb_protect_name => psb_d_par_csr_spspmm +!!$ Implicit None +!!$ +!!$ type(psb_d_csr_sparse_mat),intent(in) :: acsr +!!$ type(psb_d_csr_sparse_mat),intent(inout) :: bcsr +!!$ type(psb_d_csr_sparse_mat),intent(out) :: ccsr +!!$ type(psb_desc_type),intent(in) :: desc_a +!!$ type(psb_desc_type),intent(inout) :: desc_c +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: data +!!$ ! ...local scalars.... +!!$ integer(psb_ipk_) :: ictxt, np,me +!!$ integer(psb_ipk_) :: ncol, nnz +!!$ type(psb_d_csr_sparse_mat) :: tcsr1 +!!$ logical :: update_desc_c +!!$ integer(psb_ipk_) :: debug_level, debug_unit, err_act +!!$ character(len=20) :: name, ch_err +!!$ +!!$ if(psb_get_errstatus() /= 0) return +!!$ info=psb_success_ +!!$ name='psb_d_p_csr_spspmm' +!!$ call psb_erractionsave(err_act) +!!$ if (psb_errstatus_fatal()) then +!!$ info = psb_err_internal_error_ ; goto 9999 +!!$ end if +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ ictxt = desc_a%get_context() +!!$ +!!$ call psb_info(ictxt, me, np) +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),': Start' +!!$ +!!$ update_desc_c = desc_c%is_bld() +!!$ +!!$ ! +!!$ ! This is a bit tricky. +!!$ ! DESC_A is the descriptor of (the columns of) A, and therefore +!!$ ! of the rows of B; the columns of B, in the intended usage, span +!!$ ! a different space for which we have DESC_C. +!!$ ! We are gathering the halo rows of B to multiply by A; +!!$ ! now, the columns of B would ideally be kept in +!!$ ! global numbering, so that we can call this repeatedly to accumulate +!!$ ! the product of multiple operators, and convert to local numbering +!!$ ! at the last possible moment. However, this would imply calling +!!$ ! the serial SPSPMM with a matrix B with the GLOBAL number of columns +!!$ ! and this could be very expensive in memory. The solution is to keep B +!!$ ! in local numbering, so that only columns really appearing count, but to +!!$ ! expand the descriptor when gathering the halo, because by performing +!!$ ! the products we are extending the support of the operator; hence +!!$ ! this routine is intended to be called with a temporary descriptor +!!$ ! DESC_C which is in the BUILD state, to allow for such expansion +!!$ ! across multiple products. +!!$ ! The caller will at some later point finalize the descriptor DESC_C. +!!$ ! +!!$ +!!$ ncol = desc_a%get_local_cols() +!!$ call psb_sphalo(bcsr,desc_a,tcsr1,info,& +!!$ & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) +!!$ nnz = tcsr1%get_nzeros() +!!$ if (update_desc_c) then +!!$ call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) +!!$ else +!!$ call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) +!!$ end if +!!$ if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1) +!!$ if (info == psb_success_) call tcsr1%free() +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ call bcsr%set_ncols(desc_c%get_local_cols()) +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm 3' +!!$ if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() +!!$ call psb_spspmm(acsr,bcsr,ccsr,info) +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(ictxt,err_act) +!!$ +!!$ return +!!$ +!!$End Subroutine psb_d_par_csr_spspmm + +Subroutine psb_ld_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + use psb_base_mod, psb_protect_name => psb_ld_par_csr_spspmm + Implicit None + + type(psb_ld_csr_sparse_mat),intent(in) :: acsr + type(psb_ld_csr_sparse_mat),intent(inout) :: bcsr + type(psb_ld_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + ! ...local scalars.... + integer(psb_ipk_) :: ictxt, np,me + integer(psb_lpk_) :: nacol, nccol, nnz + type(psb_ld_csr_sparse_mat) :: tcsr1 + logical :: update_desc_c + integer(psb_ipk_) :: debug_level, debug_unit, err_act + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='psb_ld_p_csr_spspmm' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_a%get_context() + + call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + update_desc_c = desc_c%is_bld() + + ! + ! This is a bit tricky. + ! DESC_A is the descriptor of (the columns of) A, and therefore + ! of the rows of B; the columns of B, in the intended usage, span + ! a different space for which we have DESC_C. + ! We are gathering the halo rows of B to multiply by A; + ! now, the columns of B would ideally be kept in + ! global numbering, so that we can call this repeatedly to accumulate + ! the product of multiple operators, and convert to local numbering + ! at the last possible moment. However, this would imply calling + ! the serial SPSPMM with a matrix B with the GLOBAL number of columns + ! and this could be very expensive in memory. The solution is to keep B + ! in local numbering, so that only columns really appearing count, but to + ! expand the descriptor when gathering the halo, because by performing + ! the products we are extending the support of the operator; hence + ! this routine is intended to be called with a temporary descriptor + ! DESC_C which is in the BUILD state, to allow for such expansion + ! across multiple products. + ! The caller will at some later point finalize the descriptor DESC_C. + ! + + nacol = desc_a%get_local_cols() + call psb_sphalo(bcsr,desc_a,tcsr1,info,& + & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) + nnz = tcsr1%get_nzeros() + if (update_desc_c) then + call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) + else + call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) + end if + if (info == psb_success_) call psb_rwextd(nacol,bcsr,info,b=tcsr1) + if (info == psb_success_) call tcsr1%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + nccol = desc_c%get_local_cols() + call bcsr%set_ncols(nccol) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm 3' + if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() + call psb_spspmm(acsr,bcsr,ccsr,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +End Subroutine psb_ld_par_csr_spspmm diff --git a/base/tools/psb_s_par_csr_spspmm.f90 b/base/tools/psb_s_par_csr_spspmm.f90 new file mode 100644 index 00000000..b998ffc6 --- /dev/null +++ b/base/tools/psb_s_par_csr_spspmm.f90 @@ -0,0 +1,249 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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: psb_s_par_csr_spspmm.f90 +! +! Subroutine: psb_s_par_csr_spspmm +! Version: real +! +! This routine computes a parallel product of two sparse matrices +! +! C = A * B +! +! where all the matrices are stored in CSR. On input and output the matrices +! are stored with column indices in local numbering, but inermediate quantities +! are in global numbering because gathering the halo of B to multiply it +! by A implies a potential enlargement of the support. +! Also, B may have a column index space different from its row index space, +! which is obviously the same as the column space of A. +! +! +! Arguments: +! acsr - type(psb_s_csr_sparse_mat), input. +! The sparse matrix structure A +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the column space of A +! bcsr - type(psb_s_csr_sparse_mat), input/output. +! The sparse matrix structure B, gets row-extended on output +! ccsr - type(psb_s_csr_sparse_mat), output +! The sparse matrix structure C +! desc_c - type(psb_desc_type), input/output. +! The communication descriptor of the column space of B +! +! info - integer, output. +! Error code. +! +!!$Subroutine psb_s_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) +!!$ use psb_base_mod, psb_protect_name => psb_s_par_csr_spspmm +!!$ Implicit None +!!$ +!!$ type(psb_s_csr_sparse_mat),intent(in) :: acsr +!!$ type(psb_s_csr_sparse_mat),intent(inout) :: bcsr +!!$ type(psb_s_csr_sparse_mat),intent(out) :: ccsr +!!$ type(psb_desc_type),intent(in) :: desc_a +!!$ type(psb_desc_type),intent(inout) :: desc_c +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: data +!!$ ! ...local scalars.... +!!$ integer(psb_ipk_) :: ictxt, np,me +!!$ integer(psb_ipk_) :: ncol, nnz +!!$ type(psb_s_csr_sparse_mat) :: tcsr1 +!!$ logical :: update_desc_c +!!$ integer(psb_ipk_) :: debug_level, debug_unit, err_act +!!$ character(len=20) :: name, ch_err +!!$ +!!$ if(psb_get_errstatus() /= 0) return +!!$ info=psb_success_ +!!$ name='psb_s_p_csr_spspmm' +!!$ call psb_erractionsave(err_act) +!!$ if (psb_errstatus_fatal()) then +!!$ info = psb_err_internal_error_ ; goto 9999 +!!$ end if +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ ictxt = desc_a%get_context() +!!$ +!!$ call psb_info(ictxt, me, np) +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),': Start' +!!$ +!!$ update_desc_c = desc_c%is_bld() +!!$ +!!$ ! +!!$ ! This is a bit tricky. +!!$ ! DESC_A is the descriptor of (the columns of) A, and therefore +!!$ ! of the rows of B; the columns of B, in the intended usage, span +!!$ ! a different space for which we have DESC_C. +!!$ ! We are gathering the halo rows of B to multiply by A; +!!$ ! now, the columns of B would ideally be kept in +!!$ ! global numbering, so that we can call this repeatedly to accumulate +!!$ ! the product of multiple operators, and convert to local numbering +!!$ ! at the last possible moment. However, this would imply calling +!!$ ! the serial SPSPMM with a matrix B with the GLOBAL number of columns +!!$ ! and this could be very expensive in memory. The solution is to keep B +!!$ ! in local numbering, so that only columns really appearing count, but to +!!$ ! expand the descriptor when gathering the halo, because by performing +!!$ ! the products we are extending the support of the operator; hence +!!$ ! this routine is intended to be called with a temporary descriptor +!!$ ! DESC_C which is in the BUILD state, to allow for such expansion +!!$ ! across multiple products. +!!$ ! The caller will at some later point finalize the descriptor DESC_C. +!!$ ! +!!$ +!!$ ncol = desc_a%get_local_cols() +!!$ call psb_sphalo(bcsr,desc_a,tcsr1,info,& +!!$ & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) +!!$ nnz = tcsr1%get_nzeros() +!!$ if (update_desc_c) then +!!$ call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) +!!$ else +!!$ call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) +!!$ end if +!!$ if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1) +!!$ if (info == psb_success_) call tcsr1%free() +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ call bcsr%set_ncols(desc_c%get_local_cols()) +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm 3' +!!$ if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() +!!$ call psb_spspmm(acsr,bcsr,ccsr,info) +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(ictxt,err_act) +!!$ +!!$ return +!!$ +!!$End Subroutine psb_s_par_csr_spspmm + +Subroutine psb_ls_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + use psb_base_mod, psb_protect_name => psb_ls_par_csr_spspmm + Implicit None + + type(psb_ls_csr_sparse_mat),intent(in) :: acsr + type(psb_ls_csr_sparse_mat),intent(inout) :: bcsr + type(psb_ls_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + ! ...local scalars.... + integer(psb_ipk_) :: ictxt, np,me + integer(psb_lpk_) :: nacol, nccol, nnz + type(psb_ls_csr_sparse_mat) :: tcsr1 + logical :: update_desc_c + integer(psb_ipk_) :: debug_level, debug_unit, err_act + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='psb_ls_p_csr_spspmm' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_a%get_context() + + call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + update_desc_c = desc_c%is_bld() + + ! + ! This is a bit tricky. + ! DESC_A is the descriptor of (the columns of) A, and therefore + ! of the rows of B; the columns of B, in the intended usage, span + ! a different space for which we have DESC_C. + ! We are gathering the halo rows of B to multiply by A; + ! now, the columns of B would ideally be kept in + ! global numbering, so that we can call this repeatedly to accumulate + ! the product of multiple operators, and convert to local numbering + ! at the last possible moment. However, this would imply calling + ! the serial SPSPMM with a matrix B with the GLOBAL number of columns + ! and this could be very expensive in memory. The solution is to keep B + ! in local numbering, so that only columns really appearing count, but to + ! expand the descriptor when gathering the halo, because by performing + ! the products we are extending the support of the operator; hence + ! this routine is intended to be called with a temporary descriptor + ! DESC_C which is in the BUILD state, to allow for such expansion + ! across multiple products. + ! The caller will at some later point finalize the descriptor DESC_C. + ! + + nacol = desc_a%get_local_cols() + call psb_sphalo(bcsr,desc_a,tcsr1,info,& + & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) + nnz = tcsr1%get_nzeros() + if (update_desc_c) then + call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) + else + call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) + end if + if (info == psb_success_) call psb_rwextd(nacol,bcsr,info,b=tcsr1) + if (info == psb_success_) call tcsr1%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + nccol = desc_c%get_local_cols() + call bcsr%set_ncols(nccol) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm 3' + if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() + call psb_spspmm(acsr,bcsr,ccsr,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +End Subroutine psb_ls_par_csr_spspmm diff --git a/base/tools/psb_z_par_csr_spspmm.f90 b/base/tools/psb_z_par_csr_spspmm.f90 new file mode 100644 index 00000000..71e1b7ee --- /dev/null +++ b/base/tools/psb_z_par_csr_spspmm.f90 @@ -0,0 +1,249 @@ +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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: psb_z_par_csr_spspmm.f90 +! +! Subroutine: psb_z_par_csr_spspmm +! Version: complex +! +! This routine computes a parallel product of two sparse matrices +! +! C = A * B +! +! where all the matrices are stored in CSR. On input and output the matrices +! are stored with column indices in local numbering, but inermediate quantities +! are in global numbering because gathering the halo of B to multiply it +! by A implies a potential enlargement of the support. +! Also, B may have a column index space different from its row index space, +! which is obviously the same as the column space of A. +! +! +! Arguments: +! acsr - type(psb_z_csr_sparse_mat), input. +! The sparse matrix structure A +! desc_a - type(psb_desc_type), input. +! The communication descriptor of the column space of A +! bcsr - type(psb_z_csr_sparse_mat), input/output. +! The sparse matrix structure B, gets row-extended on output +! ccsr - type(psb_z_csr_sparse_mat), output +! The sparse matrix structure C +! desc_c - type(psb_desc_type), input/output. +! The communication descriptor of the column space of B +! +! info - integer, output. +! Error code. +! +!!$Subroutine psb_z_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) +!!$ use psb_base_mod, psb_protect_name => psb_z_par_csr_spspmm +!!$ Implicit None +!!$ +!!$ type(psb_z_csr_sparse_mat),intent(in) :: acsr +!!$ type(psb_z_csr_sparse_mat),intent(inout) :: bcsr +!!$ type(psb_z_csr_sparse_mat),intent(out) :: ccsr +!!$ type(psb_desc_type),intent(in) :: desc_a +!!$ type(psb_desc_type),intent(inout) :: desc_c +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_), intent(in), optional :: data +!!$ ! ...local scalars.... +!!$ integer(psb_ipk_) :: ictxt, np,me +!!$ integer(psb_ipk_) :: ncol, nnz +!!$ type(psb_z_csr_sparse_mat) :: tcsr1 +!!$ logical :: update_desc_c +!!$ integer(psb_ipk_) :: debug_level, debug_unit, err_act +!!$ character(len=20) :: name, ch_err +!!$ +!!$ if(psb_get_errstatus() /= 0) return +!!$ info=psb_success_ +!!$ name='psb_z_p_csr_spspmm' +!!$ call psb_erractionsave(err_act) +!!$ if (psb_errstatus_fatal()) then +!!$ info = psb_err_internal_error_ ; goto 9999 +!!$ end if +!!$ debug_unit = psb_get_debug_unit() +!!$ debug_level = psb_get_debug_level() +!!$ +!!$ ictxt = desc_a%get_context() +!!$ +!!$ call psb_info(ictxt, me, np) +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),': Start' +!!$ +!!$ update_desc_c = desc_c%is_bld() +!!$ +!!$ ! +!!$ ! This is a bit tricky. +!!$ ! DESC_A is the descriptor of (the columns of) A, and therefore +!!$ ! of the rows of B; the columns of B, in the intended usage, span +!!$ ! a different space for which we have DESC_C. +!!$ ! We are gathering the halo rows of B to multiply by A; +!!$ ! now, the columns of B would ideally be kept in +!!$ ! global numbering, so that we can call this repeatedly to accumulate +!!$ ! the product of multiple operators, and convert to local numbering +!!$ ! at the last possible moment. However, this would imply calling +!!$ ! the serial SPSPMM with a matrix B with the GLOBAL number of columns +!!$ ! and this could be very expensive in memory. The solution is to keep B +!!$ ! in local numbering, so that only columns really appearing count, but to +!!$ ! expand the descriptor when gathering the halo, because by performing +!!$ ! the products we are extending the support of the operator; hence +!!$ ! this routine is intended to be called with a temporary descriptor +!!$ ! DESC_C which is in the BUILD state, to allow for such expansion +!!$ ! across multiple products. +!!$ ! The caller will at some later point finalize the descriptor DESC_C. +!!$ ! +!!$ +!!$ ncol = desc_a%get_local_cols() +!!$ call psb_sphalo(bcsr,desc_a,tcsr1,info,& +!!$ & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) +!!$ nnz = tcsr1%get_nzeros() +!!$ if (update_desc_c) then +!!$ call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) +!!$ else +!!$ call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) +!!$ end if +!!$ if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1) +!!$ if (info == psb_success_) call tcsr1%free() +!!$ if(info /= psb_success_) then +!!$ call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') +!!$ goto 9999 +!!$ end if +!!$ call bcsr%set_ncols(desc_c%get_local_cols()) +!!$ +!!$ +!!$ if (debug_level >= psb_debug_outer_) & +!!$ & write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm 3' +!!$ if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& +!!$ & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() +!!$ call psb_spspmm(acsr,bcsr,ccsr,info) +!!$ +!!$ call psb_erractionrestore(err_act) +!!$ return +!!$ +!!$9999 call psb_error_handler(ictxt,err_act) +!!$ +!!$ return +!!$ +!!$End Subroutine psb_z_par_csr_spspmm + +Subroutine psb_lz_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data) + use psb_base_mod, psb_protect_name => psb_lz_par_csr_spspmm + Implicit None + + type(psb_lz_csr_sparse_mat),intent(in) :: acsr + type(psb_lz_csr_sparse_mat),intent(inout) :: bcsr + type(psb_lz_csr_sparse_mat),intent(out) :: ccsr + type(psb_desc_type),intent(in) :: desc_a + type(psb_desc_type),intent(inout) :: desc_c + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: data + ! ...local scalars.... + integer(psb_ipk_) :: ictxt, np,me + integer(psb_lpk_) :: nacol, nccol, nnz + type(psb_lz_csr_sparse_mat) :: tcsr1 + logical :: update_desc_c + integer(psb_ipk_) :: debug_level, debug_unit, err_act + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=psb_success_ + name='psb_lz_p_csr_spspmm' + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + ictxt = desc_a%get_context() + + call psb_info(ictxt, me, np) + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),': Start' + + update_desc_c = desc_c%is_bld() + + ! + ! This is a bit tricky. + ! DESC_A is the descriptor of (the columns of) A, and therefore + ! of the rows of B; the columns of B, in the intended usage, span + ! a different space for which we have DESC_C. + ! We are gathering the halo rows of B to multiply by A; + ! now, the columns of B would ideally be kept in + ! global numbering, so that we can call this repeatedly to accumulate + ! the product of multiple operators, and convert to local numbering + ! at the last possible moment. However, this would imply calling + ! the serial SPSPMM with a matrix B with the GLOBAL number of columns + ! and this could be very expensive in memory. The solution is to keep B + ! in local numbering, so that only columns really appearing count, but to + ! expand the descriptor when gathering the halo, because by performing + ! the products we are extending the support of the operator; hence + ! this routine is intended to be called with a temporary descriptor + ! DESC_C which is in the BUILD state, to allow for such expansion + ! across multiple products. + ! The caller will at some later point finalize the descriptor DESC_C. + ! + + nacol = desc_a%get_local_cols() + call psb_sphalo(bcsr,desc_a,tcsr1,info,& + & colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data) + nnz = tcsr1%get_nzeros() + if (update_desc_c) then + call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info) + else + call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info) + end if + if (info == psb_success_) call psb_rwextd(nacol,bcsr,info,b=tcsr1) + if (info == psb_success_) call tcsr1%free() + if(info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') + goto 9999 + end if + nccol = desc_c%get_local_cols() + call bcsr%set_ncols(nccol) + + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm 3' + if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),& + & 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols() + call psb_spspmm(acsr,bcsr,ccsr,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +End Subroutine psb_lz_par_csr_spspmm diff --git a/configure.ac b/configure.ac index b373eea5..18dbfdd3 100755 --- a/configure.ac +++ b/configure.ac @@ -479,11 +479,6 @@ dnl FDEFINES="$psblas_cv_define_prepend-DMPI_MOD_F08 $FDEFINES"; fi fi -dnl PAC_ARG_LONG_INTEGERS -dnl if test x"$pac_cv_long_integers" == x"yes" ; then -dnl FDEFINES="$psblas_cv_define_prepend-DLONG_INTEGERS $FDEFINES"; -dnl CDEFINES="-DLONG_INTEGERS_ $CDEFINES"; -dnl fi PAC_ARG_WITH_IPK PAC_ARG_WITH_LPK # Defaults for IPK/LPK