diff --git a/Makefile b/Makefile index d3c36e25c..ee90a48cc 100644 --- a/Makefile +++ b/Makefile @@ -98,6 +98,7 @@ clean: cleanlib cleantest: cd test/fileread && $(MAKE) clean cd test/pdegen && $(MAKE) clean + cd test/nested && $(MAKE) clean cd test/util && $(MAKE) clean cleanlib: diff --git a/base/CMakeLists.txt b/base/CMakeLists.txt index 5b01ebf76..d60799367 100644 --- a/base/CMakeLists.txt +++ b/base/CMakeLists.txt @@ -569,6 +569,7 @@ set(PSB_base_source_files modules/tools/psb_d_tools_mod.F90 modules/tools/psb_cd_nest_tools_mod.F90 modules/tools/psb_d_nest_tools_mod.F90 + modules/tools/psb_d_nest_builder_mod.F90 modules/tools/psb_c_tools_mod.F90 modules/tools/psb_e_tools_a_mod.f90 modules/tools/psb_i2_tools_a_mod.f90 @@ -641,10 +642,8 @@ set(PSB_base_source_files modules/desc/psb_glist_map_mod.F90 modules/psb_base_mod.f90 modules/desc/psb_desc_nest_mod.f90 - modules/serial/psb_d_nest_vect_mod.f90 modules/serial/psb_d_nest_mat_mod.f90 - modules/comm/psb_d_nest_comm_mod.f90 - modules/psblas/psb_d_nest_psblas_mod.f90 + modules/serial/psb_d_nest_base_mat_mod.F90 modules/psb_d_nest_mod.f90 ) foreach(file IN LISTS PSB_base_source_files) diff --git a/base/modules/Makefile b/base/modules/Makefile index bc8ecfbe8..b061d79de 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -84,7 +84,7 @@ SERIAL_MODS=serial/psb_s_serial_mod.o serial/psb_d_serial_mod.o \ serial/psb_d_base_mat_mod.o serial/psb_d_csr_mat_mod.o serial/psb_d_csc_mat_mod.o serial/psb_d_mat_mod.o \ serial/psb_c_base_mat_mod.o serial/psb_c_csr_mat_mod.o serial/psb_c_csc_mat_mod.o serial/psb_c_mat_mod.o \ serial/psb_z_base_mat_mod.o serial/psb_z_csr_mat_mod.o serial/psb_z_csc_mat_mod.o serial/psb_z_mat_mod.o \ - serial/psb_d_nest_vect_mod.o serial/psb_d_nest_mat_mod.o + serial/psb_d_nest_mat_mod.o serial/psb_d_nest_base_mat_mod.o #\ # serial/psb_ls_csr_mat_mod.o serial/psb_ld_csr_mat_mod.o serial/psb_lc_csr_mat_mod.o serial/psb_lz_csr_mat_mod.o #\ @@ -99,6 +99,7 @@ UTIL_MODS = desc/psb_desc_const_mod.o desc/psb_indx_map_mod.o\ tools/psb_i_tools_mod.o tools/psb_l_tools_mod.o \ tools/psb_s_tools_mod.o tools/psb_d_tools_mod.o\ tools/psb_d_nest_tools_mod.o \ + tools/psb_d_nest_builder_mod.o \ tools/psb_c_tools_mod.o tools/psb_z_tools_mod.o \ tools/psb_i2_tools_a_mod.o tools/psb_m_tools_a_mod.o tools/psb_e_tools_a_mod.o \ tools/psb_s_tools_a_mod.o tools/psb_d_tools_a_mod.o\ @@ -114,7 +115,6 @@ UTIL_MODS = desc/psb_desc_const_mod.o desc/psb_indx_map_mod.o\ comm/psb_s_comm_mod.o comm/psb_d_comm_mod.o\ comm/psb_c_comm_mod.o comm/psb_z_comm_mod.o \ comm/psb_i2_comm_a_mod.o \ - comm/psb_d_nest_comm_mod.o \ comm/psb_m_comm_a_mod.o comm/psb_e_comm_a_mod.o \ comm/psb_s_comm_a_mod.o comm/psb_d_comm_a_mod.o\ comm/psb_c_comm_a_mod.o comm/psb_z_comm_a_mod.o \ @@ -128,7 +128,6 @@ UTIL_MODS = desc/psb_desc_const_mod.o desc/psb_indx_map_mod.o\ psblas/psb_s_psblas_mod.o psblas/psb_c_psblas_mod.o \ psblas/psb_d_psblas_mod.o psblas/psb_z_psblas_mod.o \ psblas/psb_psblas_mod.o \ - psblas/psb_d_nest_psblas_mod.o \ psb_check_mod.o desc/psb_hash_mod.o @@ -430,9 +429,14 @@ tools/psb_cd_nest_tools_mod.o: tools/psb_cd_nest_tools_mod.F90 tools/psb_cd_tool $(FC) $(FCOPT) $(FINCLUDES) $(FDEFINES) -c tools/psb_cd_nest_tools_mod.F90 -o tools/psb_cd_nest_tools_mod.o tools/psb_d_nest_tools_mod.o: tools/psb_d_nest_tools_mod.F90 tools/psb_d_tools_mod.o \ - desc/psb_desc_nest_mod.o serial/psb_d_nest_mat_mod.o serial/psb_d_nest_vect_mod.o + desc/psb_desc_nest_mod.o serial/psb_d_nest_mat_mod.o $(FC) $(FCOPT) $(FINCLUDES) $(FDEFINES) -c tools/psb_d_nest_tools_mod.F90 -o tools/psb_d_nest_tools_mod.o +tools/psb_d_nest_builder_mod.o: tools/psb_d_nest_builder_mod.F90 tools/psb_cd_tools_mod.o \ + tools/psb_cd_nest_tools_mod.o tools/psb_d_nest_tools_mod.o \ + serial/psb_d_nest_base_mat_mod.o serial/psb_d_nest_mat_mod.o desc/psb_desc_nest_mod.o + $(FC) $(FCOPT) $(FINCLUDES) $(FDEFINES) -c tools/psb_d_nest_builder_mod.F90 -o tools/psb_d_nest_builder_mod.o + tools/psb_cd_tools_mod.o tools/psb_i_tools_mod.o tools/psb_l_tools_mod.o \ tools/psb_s_tools_mod.o tools/psb_d_tools_mod.o \ tools/psb_c_tools_mod.o tools/psb_z_tools_mod.o \ @@ -454,26 +458,17 @@ psblas/psb_z_psblas_mod.o: serial/psb_z_vect_mod.o serial/psb_z_mat_mod.o psblas/psb_psblas_mod.o: psblas/psb_s_psblas_mod.o psblas/psb_c_psblas_mod.o psblas/psb_d_psblas_mod.o psblas/psb_z_psblas_mod.o psblas/psb_s_psblas_mod.o psblas/psb_c_psblas_mod.o psblas/psb_d_psblas_mod.o psblas/psb_z_psblas_mod.o: serial/psb_mat_mod.o desc/psb_desc_mod.o -# --- nested mat/vec/desc dependencies --- +# --- nested mat/desc dependencies (MATNEST) --- desc/psb_desc_nest_mod.o: desc/psb_desc_mod.o -serial/psb_d_nest_vect_mod.o: serial/psb_d_vect_mod.o desc/psb_desc_mod.o serial/psb_d_nest_mat_mod.o: serial/psb_d_mat_mod.o -comm/psb_d_nest_comm_mod.o: \ - desc/psb_desc_nest_mod.o \ - serial/psb_d_nest_vect_mod.o \ - comm/psb_d_comm_mod.o -psblas/psb_d_nest_psblas_mod.o: \ - desc/psb_desc_nest_mod.o \ - serial/psb_d_nest_vect_mod.o \ - serial/psb_d_nest_mat_mod.o \ - serial/psb_d_mat_mod.o \ - psblas/psb_d_psblas_mod.o \ - comm/psb_d_nest_comm_mod.o +serial/psb_d_nest_base_mat_mod.o: serial/psb_d_nest_mat_mod.o desc/psb_desc_nest_mod.o serial/psb_d_base_mat_mod.o serial/psb_d_mat_mod.o desc/psb_desc_mod.o psb_d_nest_mod.o: \ desc/psb_desc_nest_mod.o \ - serial/psb_d_nest_vect_mod.o \ serial/psb_d_nest_mat_mod.o \ - psblas/psb_d_nest_psblas_mod.o + serial/psb_d_nest_base_mat_mod.o \ + tools/psb_cd_nest_tools_mod.o \ + tools/psb_d_nest_tools_mod.o \ + tools/psb_d_nest_builder_mod.o psb_base_mod.o: $(MODULES) diff --git a/base/modules/comm/psb_d_nest_comm_mod.f90 b/base/modules/comm/psb_d_nest_comm_mod.f90 deleted file mode 100644 index cef64ebf7..000000000 --- a/base/modules/comm/psb_d_nest_comm_mod.f90 +++ /dev/null @@ -1,93 +0,0 @@ -! -! 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. -! -! -! module: psb_d_nest_comm_mod -! -! Communication operations for nested (block-structured) double precision -! real vectors. -! -! psb_d_nest_halo -! Halo exchange for all column blocks of a nested vector. -! Calls psb_halo(x(j), descs(1,j)) for each column block j. -! All descriptors descs(i,j) for fixed j are equivalent; -! Called once before block SpMM to populate ghost entries of x. -! -! psb_d_nest_ovrl -! Overlap update for all row blocks of a nested vector. -! Calls psb_ovrl(x(i), descs(i,i)) for each row block i using the -! diagonal descriptor. -! Called after operations that contribute to overlapping rows -! (e.g. FEM assembly). -! -module psb_d_nest_comm_mod - use psb_desc_nest_mod - use psb_d_nest_vect_mod - use psb_d_comm_mod, only : psb_halo, psb_ovrl - use psb_const_mod, only : psb_ipk_ - implicit none - - private - public :: psb_d_nest_halo, psb_d_nest_ovrl - -contains - - subroutine psb_d_nest_halo(xnest, descs, info, tran, mode, data) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - character, optional, intent(in) :: tran - integer(psb_ipk_), optional, intent(in) :: mode, data - - integer(psb_ipk_) :: j - - info = 0 - do j = 1, xnest%nblocks - call psb_halo(xnest%vects(j), descs%descs(1,j), info, tran=tran, mode=mode, data=data) - if (info /= 0) return - end do - end subroutine psb_d_nest_halo - - subroutine psb_d_nest_ovrl(xnest, descs, info, update, mode) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_), optional, intent(in) :: update, mode - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_ovrl(xnest%vects(i), descs%descs(i,i), info, update=update, mode=mode) - if (info /= 0) return - end do - end subroutine psb_d_nest_ovrl - -end module psb_d_nest_comm_mod diff --git a/base/modules/desc/psb_desc_nest_mod.f90 b/base/modules/desc/psb_desc_nest_mod.f90 index 5260285bf..d74d27869 100644 --- a/base/modules/desc/psb_desc_nest_mod.f90 +++ b/base/modules/desc/psb_desc_nest_mod.f90 @@ -30,6 +30,7 @@ ! ! ! module: psb_desc_nest_mod +! Author: Simone Staccone (Stack-1) ! ! Defines psb_desc_nest_type: a 2-D array of psb_desc_type objects, ! one per block matrix entry in an nrblocks x ncblocks block system. @@ -37,6 +38,8 @@ ! module psb_desc_nest_mod use psb_desc_mod + use psb_error_mod + implicit none type :: psb_desc_nest_type @@ -69,31 +72,48 @@ contains end function psb_desc_nest_get_ncblocks ! get_desc: copy descriptor (i,j) into the output argument - subroutine psb_desc_nest_get_desc(d, i, j, desc, info) + subroutine psb_desc_nest_get_desc(d, i_block_row, j_block_col, desc, info) class(psb_desc_nest_type), intent(in) :: d - integer(psb_ipk_), intent(in) :: i, j + integer(psb_ipk_), intent(in) :: i_block_row, j_block_col type(psb_desc_type), intent(out):: desc integer(psb_ipk_), intent(out):: info + character(len=64) :: name info = 0 - if (i < 1 .or. i > d%nrblocks .or. j < 1 .or. j > d%ncblocks) then + name = 'psb_desc_nest_get_desc' + + if (i_block_row < 1 .or. i_block_row > d%nrblocks .or. & + & j_block_col < 1 .or. j_block_col > d%ncblocks) then info = -1 + call psb_errpush(info, name, a_err='Invalid block indices') return end if - desc = d%descs(i,j) + desc = d%descs(i_block_row,j_block_col) end subroutine psb_desc_nest_get_desc - ! is_valid: true if all diagonal sub-descriptors are valid + ! is_valid: true if the per-column descriptors used by the kernel are valid. + ! The previous version only checked the diagonal descs(i,i), which is + ! wrong for saddle-point systems where the (2,2) block (and hence its + ! diagonal descriptor) is absent. The nested halo relies on the per-column + ! descriptors descs(1,j) (all descs(i,j) for fixed j are equivalent), so we + ! validate those instead of the diagonal. function psb_desc_nest_is_valid(d) result(valid) class(psb_desc_nest_type), intent(in) :: d - logical :: valid - integer(psb_ipk_) :: i + logical :: valid + integer(psb_ipk_) :: j_block_col, info + character(len=64) :: name + character(len=20) :: colid + + name = 'psb_desc_nest_is_valid' + info = 0 valid = (d%nrblocks >= 1) .and. (d%ncblocks >= 1) .and. allocated(d%descs) if (valid) then - do i = 1, min(d%nrblocks, d%ncblocks) - if (.not. d%descs(i,i)%is_valid()) then + do j_block_col = 1, d%ncblocks + if (.not. d%descs(1,j_block_col)%is_valid()) then valid = .false. + info = -1 + call psb_errpush(info, name, a_err='Invalid descriptor in column '//trim(colid)) return end if end do @@ -101,16 +121,16 @@ contains end function psb_desc_nest_is_valid ! sizeof: total memory (bytes) of all sub-descriptors - function psb_desc_nest_sizeof(d) result(s) + function psb_desc_nest_sizeof(d) result(total_bytes) class(psb_desc_nest_type), intent(in) :: d - integer(psb_epk_) :: s - integer(psb_ipk_) :: i, j + integer(psb_epk_) :: total_bytes + integer(psb_ipk_) :: i_block_row, j_block_col - s = 0_psb_epk_ + total_bytes = 0_psb_epk_ if (allocated(d%descs)) then - do j = 1, d%ncblocks - do i = 1, d%nrblocks - s = s + d%descs(i,j)%sizeof() + do j_block_col = 1, d%ncblocks + do i_block_row = 1, d%nrblocks + total_bytes = total_bytes + d%descs(i_block_row,j_block_col)%sizeof() end do end do end if @@ -121,18 +141,18 @@ contains class(psb_desc_nest_type), intent(inout) :: d integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, j, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, local_info info = 0 if (allocated(d%descs)) then - do j = 1, d%ncblocks - do i = 1, d%nrblocks - call d%descs(i,j)%free(linfo) - if (linfo /= 0 .and. info == 0) info = linfo + do j_block_col = 1, d%ncblocks + do i_block_row = 1, d%nrblocks + call d%descs(i_block_row,j_block_col)%free(local_info) + if (local_info /= 0 .and. info == 0) info = local_info end do end do - deallocate(d%descs, stat=linfo) - if (linfo /= 0 .and. info == 0) info = linfo + deallocate(d%descs, stat=local_info) + if (local_info /= 0 .and. info == 0) info = local_info end if d%nrblocks = 0 d%ncblocks = 0 diff --git a/base/modules/psb_d_nest_mod.f90 b/base/modules/psb_d_nest_mod.f90 index 93f4432c1..367055100 100644 --- a/base/modules/psb_d_nest_mod.f90 +++ b/base/modules/psb_d_nest_mod.f90 @@ -30,6 +30,7 @@ ! ! ! module: psb_d_nest_mod +! Author: Simone Staccone (Stack-1) ! ! Umbrella module for the nested (block-structured) double precision ! real types. Users need only: @@ -39,8 +40,10 @@ ! to access all three container types and their parallel operations. ! module psb_d_nest_mod - use psb_desc_nest_mod - use psb_d_nest_vect_mod - use psb_d_nest_mat_mod - use psb_d_nest_psblas_mod + use psb_desc_nest_mod ! grid descriptor (per-field, input to compose) + use psb_d_nest_mat_mod ! block storage (psb_d_nest_sparse_mat) + use psb_d_nest_base_mat_mod ! MATNEST operator + field-split interface + use psb_cd_nest_tools_mod ! psb_cd_nest_compose (global descriptor) + use psb_d_nest_tools_mod ! block assembly + psb_d_nest_rect_block + use psb_d_nest_builder_mod ! psb_d_nest_matrix: init/ins/asb frontend end module psb_d_nest_mod diff --git a/base/modules/psblas/psb_d_nest_psblas_mod.f90 b/base/modules/psblas/psb_d_nest_psblas_mod.f90 deleted file mode 100644 index 240065bdc..000000000 --- a/base/modules/psblas/psb_d_nest_psblas_mod.f90 +++ /dev/null @@ -1,626 +0,0 @@ -! -! 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. -! -! -! module: psb_d_nest_psblas_mod -! -! Parallel BLAS operations for the nested (block-structured) double -! precision real types. -! -! psb_d_nest_spmm -! Computes y = alpha * A_nest * x + beta * y (block SpMV). -! Three-phase algorithm: -! Phase 1 — scale y upfront: -! if beta == 0: zero all y(i) -! elif beta /= 1: y(i) = beta * y(i) for each block i -! Phase 2 — single halo exchange per column block: -! call psb_d_nest_halo(xnest, descs) -! Populates ghost entries of x(j) using descs(1,j) for each j. -! All descs(i,j) for fixed j share the same column space, so -! one exchange covers all block-rows. -! Phase 3 — local SpMM accumulation (no further communication): -! For each present block (i,j): -! y(i) += alpha * A(i,j) * x(j) -! (psb_spmm called with doswap=.false. to skip internal halo) -! -! psb_d_nest_geaxpby -! Computes y(i) = alpha * x(i) + beta * y(i) for each block i. -! -! psb_d_nest_genrm2 -! Computes ||x||_2 = sqrt( sum_i ||x(i)||_2^2 ) with a single -! global reduction. -! -! psb_d_nest_genrm2s -! Subroutine form of psb_d_nest_genrm2 (result via intent(out) argument). -! -! psb_d_nest_gedot -! Computes dot(x,y) = sum_i dot(x(i), y(i)) with a single global -! reduction. -! -! psb_d_nest_geamax -! Computes ||x||_inf = max_i ||x(i)||_inf with a single global reduction. -! -! psb_d_nest_geasum -! Computes ||x||_1 = sum_i ||x(i)||_1 with a single global reduction. -! -! psb_d_nest_gemin -! Computes min(x) = min_i min(x(i)) with a single global reduction. -! -! psb_d_nest_minquotient -! Computes min(x/y) = min_i min(x(i)/y(i)) with a single global reduction. -! -! psb_d_nest_gemlt -! Computes y(i) = x(i) .* y(i) element-wise for each block i. -! -! psb_d_nest_gediv -! Computes y(i) = x(i) ./ y(i) element-wise for each block i. -! -! psb_d_nest_geinv -! Computes y(i) = 1 / x(i) element-wise for each block i. -! -! psb_d_nest_geabs -! Computes y(i) = |x(i)| element-wise for each block i. -! -! psb_d_nest_geaddconst -! Computes z(i) = x(i) + b for each block i (b is a scalar). -! -! psb_d_nest_gecmp -! Computes z(i) = cmp(x(i), c) for each block i (c is a scalar). -! -! psb_d_nest_mask -! Applies mask operation to each block i; t is .true. iff all blocks -! return .true. -! -! psb_d_nest_upd_xyz -! Applies psb_upd_xyz(alpha,beta,gamma,delta, x(i),y(i),z(i)) for each block i. -! -! psb_d_nest_spsm -! Block-diagonal triangular solve: applies psb_spsm to each diagonal -! block (i,i) of tnest independently. -! -module psb_d_nest_psblas_mod - use psb_desc_nest_mod - use psb_d_nest_vect_mod - use psb_d_nest_mat_mod - use psb_d_mat_mod, only : psb_dspmat_type, psb_csmm - use psb_d_psblas_mod, only : psb_spmm, psb_geaxpby, psb_genrm2, psb_gedot, & - & psb_geamax, psb_geasum, psb_gemin, psb_minquotient, & - & psb_gemlt, psb_gediv, psb_geinv, psb_geabs, psb_geaddconst, & - & psb_gecmp, psb_mask, psb_upd_xyz, psb_spsm - use psb_d_nest_comm_mod, only : psb_d_nest_halo, psb_d_nest_ovrl - use psb_penv_mod, only : psb_sum, psb_max, psb_min, psb_info - use psb_const_mod, only : psb_dpk_, psb_ipk_, psb_epk_, psb_ctxt_type, dzero, done - implicit none - - private - public :: psb_d_nest_spmm, psb_d_nest_geaxpby, & - psb_d_nest_genrm2, psb_d_nest_genrm2s, psb_d_nest_gedot, & - psb_d_nest_geamax, psb_d_nest_geasum, psb_d_nest_gemin, & - psb_d_nest_minquotient, & - psb_d_nest_gemlt, psb_d_nest_gediv, psb_d_nest_geinv, & - psb_d_nest_halo, psb_d_nest_ovrl, & - psb_d_nest_geabs, psb_d_nest_geaddconst, psb_d_nest_gecmp, & - psb_d_nest_mask, psb_d_nest_upd_xyz, psb_d_nest_spsm - -contains - - ! y = alpha * A_nest * x + beta * y - subroutine psb_d_nest_spmm(alpha, anest, xnest, beta, ynest, descs, info, trans) - - real(psb_dpk_), intent(in) :: alpha, beta - type(psb_d_nest_sparse_mat), intent(in) :: anest - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - character, optional, intent(in) :: trans - - integer(psb_ipk_) :: i, j - character :: trans_ - - info = 0 - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - - if (beta == dzero) then - do i = 1, anest%nrblocks - call ynest%vects(i)%zero() - end do - else if (beta /= done) then - do i = 1, anest%nrblocks - call ynest%vects(i)%scal(beta) - end do - end if - - call psb_d_nest_halo(xnest, descs, info) - if (info /= 0) return - - do i = 1, anest%nrblocks - do j = 1, anest%ncblocks - if (anest%has_block(i, j)) then - ! y(i) += alpha * A(i,j) * x(j) (doswap=.false. skips internal halo, already done in Phase 2) - call psb_spmm(alpha, anest%mats(i,j), xnest%vects(j), & - & done, ynest%vects(i), descs%descs(i,j), info, trans=trans_, doswap=.false.) - if (info /= 0) return - end if - end do - end do - end subroutine psb_d_nest_spmm - - ! y(i) = alpha * x(i) + beta * y(i) for each block i - subroutine psb_d_nest_geaxpby(alpha, xnest, beta, ynest, descs, info) - real(psb_dpk_), intent(in) :: alpha, beta - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_geaxpby(alpha, xnest%vects(i), beta, ynest%vects(i), & - & descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_geaxpby - - ! ||x||_2 = sqrt( sum_i ||x(i)||_2^2 ) - ! Uses a single global MPI_Allreduce across all blocks. - ! global (optional, default .true.): if .false., skips MPI_Allreduce and returns - ! the process-local partial norm; use when the caller manages the reduction itself. - - function psb_d_nest_genrm2(xnest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - real(psb_dpk_) :: loc_sum, blk_nrm - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - loc_sum = dzero - do i = 1, xnest%nblocks - ! global=.false. returns local partial norm (sqrt of local partial sum) - blk_nrm = psb_genrm2(xnest%vects(i), descs%descs(i,i), info, global=.false.) - if (info /= 0) then - res = dzero - return - end if - loc_sum = loc_sum + blk_nrm * blk_nrm - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_sum(ctxt, loc_sum) - end if - - res = sqrt(loc_sum) - end function psb_d_nest_genrm2 - - ! psb_d_nest_gedot - ! dot(x, y) = sum_i dot(x(i), y(i)) - ! Uses a single global MPI_Allreduce across all blocks. - function psb_d_nest_gedot(xnest, ynest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - real(psb_dpk_) :: loc_sum, blk_dot - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - loc_sum = dzero - do i = 1, xnest%nblocks - blk_dot = psb_gedot(xnest%vects(i), ynest%vects(i), descs%descs(i,i), & - & info, global=.false.) - if (info /= 0) then - res = dzero - return - end if - loc_sum = loc_sum + blk_dot - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_sum(ctxt, loc_sum) - end if - - res = loc_sum - end function psb_d_nest_gedot - - ! Subroutine form: res = ||x||_2 (single global reduction) - subroutine psb_d_nest_genrm2s(res, xnest, descs, info, global) - real(psb_dpk_), intent(out) :: res - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - - res = psb_d_nest_genrm2(xnest, descs, info, global) - end subroutine psb_d_nest_genrm2s - - ! ||x||_inf = max_i ||x(i)||_inf (single global reduction) - function psb_d_nest_geamax(xnest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - real(psb_dpk_) :: blk_val - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - res = dzero - do i = 1, xnest%nblocks - blk_val = psb_geamax(xnest%vects(i), descs%descs(i,i), info, global=.false.) - if (info /= 0) return - if (blk_val > res) res = blk_val - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_max(ctxt, res) - end if - end function psb_d_nest_geamax - - ! ||x||_1 = sum_i ||x(i)||_1 (single global reduction) - ! function psb_d_nest_geasum(xnest, descs, info, global) result(res) - ! type(psb_d_nest_vect_type), intent(inout) :: xnest - ! type(psb_desc_nest_type), intent(in) :: descs - ! integer(psb_ipk_), intent(out) :: info - ! logical, optional, intent(in) :: global - ! real(psb_dpk_) :: res - - ! integer(psb_ipk_) :: i - ! real(psb_dpk_) :: blk_val - ! logical :: global_ - ! type(psb_ctxt_type) :: ctxt - - ! global_ = .true. - ! if (present(global)) global_ = global - - ! info = 0 - ! res = dzero - ! do i = 1, xnest%nblocks - ! blk_val = psb_geasum(xnest%vects(i), descs%descs(i,i), info, global=.false.) - ! if (info /= 0) return - ! res = res + blk_val - ! end do - - ! if (global_) then - ! ctxt = descs%descs(1,1)%get_context() - ! call psb_sum(ctxt, res) - ! end if - ! end function psb_d_nest_geasum - - ! ||x||_1 = sum_i ||x(i)||_1 - function psb_d_nest_geasum(xnest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - integer(psb_ipk_) :: nloc - real(psb_dpk_) :: blk_val - real(psb_dpk_), allocatable :: blk_vals(:) - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - res = dzero - do i = 1, xnest%nblocks - nloc = descs%descs(i,i)%get_local_rows() - blk_vals = xnest%vects(i)%get_vect(nloc) - if (size(blk_vals) > 0) then - blk_val = sum(abs(blk_vals)) - else - blk_val = dzero - end if - res = res + blk_val - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_sum(ctxt, res) - end if - end function psb_d_nest_geasum - - ! min(x) = min_i min(x(i)) - function psb_d_nest_gemin(xnest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - real(psb_dpk_) :: blk_val - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - res = huge(dzero) - do i = 1, xnest%nblocks - blk_val = psb_gemin(xnest%vects(i), descs%descs(i,i), info, global=.false.) - if (info /= 0) return - if (blk_val < res) res = blk_val - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_min(ctxt, res) - end if - end function psb_d_nest_gemin - - ! min(x/y) = min_i min(x(i)/y(i)) (single global reduction) - function psb_d_nest_minquotient(xnest, ynest, descs, info, global) result(res) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - logical, optional, intent(in) :: global - real(psb_dpk_) :: res - - integer(psb_ipk_) :: i - real(psb_dpk_) :: blk_val - logical :: global_ - type(psb_ctxt_type) :: ctxt - - global_ = .true. - if (present(global)) global_ = global - - info = 0 - res = huge(dzero) - do i = 1, xnest%nblocks - blk_val = psb_minquotient(xnest%vects(i), ynest%vects(i), & - & descs%descs(i,i), info, global=.false.) - if (info /= 0) return - if (blk_val < res) res = blk_val - end do - - if (global_) then - ctxt = descs%descs(1,1)%get_context() - call psb_min(ctxt, res) - end if - end function psb_d_nest_minquotient - - ! y(i) = x(i) .* y(i) element-wise for each block i - subroutine psb_d_nest_gemlt(xnest, ynest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_gemlt(xnest%vects(i), ynest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_gemlt - - ! y(i) = x(i) ./ y(i) element-wise for each block i - subroutine psb_d_nest_gediv(xnest, ynest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_gediv(xnest%vects(i), ynest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_gediv - - ! y(i) = 1/x(i) element-wise for each block i - subroutine psb_d_nest_geinv(xnest, ynest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_geinv(xnest%vects(i), ynest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_geinv - - ! y(i) = |x(i)| element-wise for each block i - subroutine psb_d_nest_geabs(xnest, ynest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_geabs(xnest%vects(i), ynest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_geabs - - ! z(i) = x(i) + b for each block i (b is a scalar) - subroutine psb_d_nest_geaddconst(xnest, b, znest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - real(psb_dpk_), intent(in) :: b - type(psb_d_nest_vect_type), intent(inout) :: znest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_geaddconst(xnest%vects(i), b, znest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_geaddconst - - ! z(i) = cmp(x(i), c) for each block i (c is a scalar) - subroutine psb_d_nest_gecmp(xnest, c, znest, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: xnest - real(psb_dpk_), intent(in) :: c - type(psb_d_nest_vect_type), intent(inout) :: znest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_gecmp(xnest%vects(i), c, znest%vects(i), descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_gecmp - - subroutine psb_d_nest_mask(cnest, xnest, mnest, t, descs, info) - type(psb_d_nest_vect_type), intent(inout) :: cnest - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: mnest - logical, intent(out) :: t - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - logical :: t_dummy - real(psb_dpk_) :: mmax - - info = 0 - do i = 1, cnest%nblocks - ! psb_mask(c, x, m, t, desc) semantics after the double-swap in - ! d_vect_mask_v → d_base_mask_v → d_base_mask_a: - ! first arg → constraint-type selector in d_base_mask_a - ! second arg → value to test in d_base_mask_a - ! So pass xnest (constraint types) first, cnest (values) second. - call psb_mask(xnest%vects(i), cnest%vects(i), mnest%vects(i), & - & t_dummy, descs%descs(i,i), info) - if (info /= 0) return - end do - ! t = .true. iff no block has any violated entry (mnest=1). - mmax = psb_d_nest_geamax(mnest, descs, info) - t = (mmax < 0.5_psb_dpk_) - end subroutine psb_d_nest_mask - - - ! Applies psb_upd_xyz(alpha,beta,gamma,delta,x(i),y(i),z(i)) - ! for each block i. - subroutine psb_d_nest_upd_xyz(alpha, beta, gamma, delta, xnest, ynest, znest, descs, info) - real(psb_dpk_), intent(in) :: alpha, beta, gamma, delta - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_d_nest_vect_type), intent(inout) :: znest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, xnest%nblocks - call psb_upd_xyz(alpha, beta, gamma, delta, & - & xnest%vects(i), ynest%vects(i), znest%vects(i), & - & descs%descs(i,i), info) - if (info /= 0) return - end do - end subroutine psb_d_nest_upd_xyz - - ! Block-diagonal triangular solve: applies psb_spsm to each - ! diagonal block (i,i) of tnest independently. - ! y(i) = alpha * T(i,i)^{-1} x(i) + beta * y(i) - subroutine psb_d_nest_spsm(alpha, tnest, xnest, beta, ynest, descs, info, & - & trans, scale, choice) - real(psb_dpk_), intent(in) :: alpha, beta - type(psb_d_nest_sparse_mat), intent(inout) :: tnest - type(psb_d_nest_vect_type), intent(inout) :: xnest - type(psb_d_nest_vect_type), intent(inout) :: ynest - type(psb_desc_nest_type), intent(in) :: descs - integer(psb_ipk_), intent(out) :: info - character, optional, intent(in) :: trans, scale - integer(psb_ipk_), optional, intent(in) :: choice - - integer(psb_ipk_) :: i - - info = 0 - do i = 1, tnest%nrblocks - if (.not. tnest%has_block(i, i)) then - ! No diagonal block: treat as identity => y(i) = alpha*x(i) + beta*y(i) - call psb_geaxpby(alpha, xnest%vects(i), beta, ynest%vects(i), & - & descs%descs(i,i), info) - else - call psb_spsm(alpha, tnest%mats(i,i), xnest%vects(i), beta, ynest%vects(i), & - & descs%descs(i,i), info, trans=trans, scale=scale, & - & choice=choice) - end if - if (info /= 0) return - end do - end subroutine psb_d_nest_spsm - -end module psb_d_nest_psblas_mod diff --git a/base/modules/serial/psb_d_nest_base_mat_mod.F90 b/base/modules/serial/psb_d_nest_base_mat_mod.F90 new file mode 100644 index 000000000..af135dcc2 --- /dev/null +++ b/base/modules/serial/psb_d_nest_base_mat_mod.F90 @@ -0,0 +1,405 @@ +! +! 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 prior 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_nest_base_mat_mod.F90 +! +! Module: psb_d_nest_base_mat_mod +! Author: Simone Staccone (Stack-1) +! +! Adapter that makes a block-structured (nested) operator look like a standard +! local sparse matrix to PSBLAS: psb_d_nest_base_mat EXTENDS +! psb_d_base_sparse_mat and implements csmv (the local matrix-vector product). +! Wrapped in a psb_dspmat_type and paired with the composed global descriptor +! (see psb_cd_nest_compose), the nested operator can then be fed to psb_spmm, +! psb_krylov and the AMG4PSBLAS preconditioners unchanged (MATNEST-style). +! +! The local vector handed to csmv lives in the GLOBAL local layout produced by +! psb_cd_nest_compose: the owned entries of all fields are concatenated, followed +! by the global halo. For each field we precompute field_map(field)%global_local_pos, +! the positions in that global local vector of the field's own local vector +! (owned entries first, then the field's ghosts), so we can gather the field +! input sub-vector and scatter the field output sub-vector without further +! communication (the halo exchange is done once by psb_spmm on the global desc). +! +module psb_d_nest_base_mat_mod + use psb_const_mod + use psb_error_mod + use psb_d_base_mat_mod, only : psb_d_base_sparse_mat + use psb_desc_mod, only : psb_desc_type + use psb_desc_nest_mod, only : psb_desc_nest_type + use psb_d_nest_mat_mod, only : psb_d_nest_sparse_mat + use psb_d_mat_mod, only : psb_dspmat_type + implicit none + + ! Per-field gather/scatter map into the global local vector. + ! global_local_pos(1 : n_owned) -> the field's owned entries + ! global_local_pos(n_owned+1 : size) -> the field's ghost (halo) entries + type :: psb_d_nest_field_map + integer(psb_ipk_) :: n_owned = 0 + integer(psb_ipk_), allocatable :: global_local_pos(:) + end type psb_d_nest_field_map + + type, extends(psb_d_base_sparse_mat) :: psb_d_nest_base_mat + integer(psb_ipk_) :: n_fields = 0 + type(psb_d_nest_sparse_mat), pointer :: block_storage => null() ! blocks (not owned) + type(psb_desc_nest_type), pointer :: grid_desc => null() ! per-field descriptors (not owned) + type(psb_d_nest_field_map), allocatable :: field_map(:) + contains + procedure, pass(a) :: csmv => psb_d_nest_base_csmv + procedure, pass(a) :: get_nzeros => psb_d_nest_base_get_nzeros + procedure, nopass :: get_fmt => psb_d_nest_base_get_fmt + procedure, pass(a) :: free => psb_d_nest_base_free + end type psb_d_nest_base_mat + + private + public :: psb_d_nest_base_mat, psb_d_nest_base_setup, psb_d_nest_apply_block + ! field-split interface (for the block preconditioner) + public :: psb_d_nest_get_n_fields, psb_d_nest_get_field_owned, & + & psb_d_nest_get_block, psb_d_nest_get_field_desc, & + & psb_d_nest_restrict_field, psb_d_nest_prolong_field + +contains + + function psb_d_nest_base_get_fmt() result(format_name) + character(len=5) :: format_name + format_name = 'NEST' + end function psb_d_nest_base_get_fmt + + ! free: the nested operator does NOT own block_storage / grid_desc (they are + ! pointers into the caller), so we only detach them and release the field maps. + subroutine psb_d_nest_base_free(a) + class(psb_d_nest_base_mat), intent(inout) :: a + a%block_storage => null() + a%grid_desc => null() + if (allocated(a%field_map)) deallocate(a%field_map) + a%n_fields = 0 + call a%set_null() + end subroutine psb_d_nest_base_free + + function psb_d_nest_base_get_nzeros(a) result(total_nzeros) + class(psb_d_nest_base_mat), intent(in) :: a + integer(psb_ipk_) :: total_nzeros + integer(psb_ipk_) :: i_block_row, j_block_col + total_nzeros = 0 + if (associated(a%block_storage)) then + do j_block_col = 1, a%block_storage%ncblocks + do i_block_row = 1, a%block_storage%nrblocks + if (a%block_storage%has_block(i_block_row, j_block_col)) & + & total_nzeros = total_nzeros + & + & a%block_storage%mats(i_block_row, j_block_col)%get_nzeros() + end do + end do + end if + end function psb_d_nest_base_get_nzeros + + ! Build the per-field gather maps and set the local dimensions, from the nested + ! grid descriptor (per-field distribution desc_grid%descs(1,field)) and the + ! composed global descriptor desc_global (produced by psb_cd_nest_compose). + subroutine psb_d_nest_base_setup(nest_op, block_storage, desc_grid, desc_global, info) + type(psb_d_nest_base_mat), intent(inout) :: nest_op + type(psb_d_nest_sparse_mat), target, intent(in) :: block_storage + type(psb_desc_nest_type), target, intent(in) :: desc_grid + type(psb_desc_type), intent(in) :: desc_global + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_fields, i_field, i_entry + integer(psb_ipk_) :: n_owned, n_local, n_ghost, owned_offset, local_pos + integer(psb_lpk_) :: global_idx + integer(psb_lpk_), allocatable :: field_global_offset(:) + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_base_setup' + + if (desc_grid%nrblocks /= desc_grid%ncblocks) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='nested block structure must be square') + return + end if + + n_fields = desc_grid%ncblocks + nest_op%n_fields = n_fields + nest_op%grid_desc => desc_grid + nest_op%block_storage => block_storage + + ! global field offsets (used to form ghost global indices) + allocate(field_global_offset(n_fields+1), nest_op%field_map(n_fields), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + field_global_offset(1) = 0 + do i_field = 1, n_fields + field_global_offset(i_field+1) = field_global_offset(i_field) & + & + desc_grid%descs(1,i_field)%get_global_rows() + end do + + owned_offset = 0 ! running owned-local offset in the global local vector + do i_field = 1, n_fields + n_owned = desc_grid%descs(1,i_field)%get_local_rows() + n_local = desc_grid%descs(1,i_field)%get_local_cols() + n_ghost = n_local - n_owned + nest_op%field_map(i_field)%n_owned = n_owned + allocate(nest_op%field_map(i_field)%global_local_pos(n_local), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + ! owned entries: contiguous in the global local vector + do i_entry = 1, n_owned + nest_op%field_map(i_field)%global_local_pos(i_entry) = owned_offset + i_entry + end do + ! ghost entries: locate the field's ghost global index in the global descriptor + do i_entry = 1, n_ghost + call desc_grid%descs(1,i_field)%l2g(n_owned + i_entry, global_idx, info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='l2g'); return + end if + call desc_global%g2l(field_global_offset(i_field) + global_idx, local_pos, info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='g2l'); return + end if + nest_op%field_map(i_field)%global_local_pos(n_owned + i_entry) = local_pos + end do + owned_offset = owned_offset + n_owned + end do + + call nest_op%set_nrows(desc_global%get_local_rows()) + call nest_op%set_ncols(desc_global%get_local_cols()) + call nest_op%set_asb() + + end subroutine psb_d_nest_base_setup + + ! Local block matrix-vector product: y = alpha * A_nest * x + beta * y. + ! x is in the global local layout (owned fields concatenated + global halo); + ! y holds the owned entries (global local rows). + subroutine psb_d_nest_base_csmv(alpha, a, x, beta, y, info, trans) + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + class(psb_d_nest_base_mat), intent(in) :: a + real(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + real(psb_dpk_), allocatable :: x_field(:), y_field(:) + integer(psb_ipk_) :: i_block_row, j_block_col, i_entry + integer(psb_ipk_) :: n_local_col_field, n_owned_row_field + character :: trans_op + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_base_csmv' + trans_op = 'N' + if (present(trans)) trans_op = trans + if (trans_op /= 'N' .and. trans_op /= 'n') then + ! Transposed nested product is not implemented (would swap block indices + ! and need the transposed halo); reject explicitly. See P5. + info = psb_err_transpose_not_n_unsupported_ + call psb_errpush(info, name) + return + end if + if (.not. associated(a%block_storage)) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='nested operator not set up') + return + end if + + ! y <- beta * y + if (beta == dzero) then + y(:) = dzero + else if (beta /= done) then + y(:) = beta * y(:) + end if + + do j_block_col = 1, a%n_fields + n_local_col_field = size(a%field_map(j_block_col)%global_local_pos) + allocate(x_field(n_local_col_field), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + ! gather the column-field input sub-vector (owned + that field's ghosts) + do i_entry = 1, n_local_col_field + x_field(i_entry) = x(a%field_map(j_block_col)%global_local_pos(i_entry)) + end do + + do i_block_row = 1, a%n_fields + if (a%block_storage%has_block(i_block_row, j_block_col)) then + n_owned_row_field = a%field_map(i_block_row)%n_owned + allocate(y_field(n_owned_row_field), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + ! current row-field output sub-vector (owned) + do i_entry = 1, n_owned_row_field + y_field(i_entry) = y(a%field_map(i_block_row)%global_local_pos(i_entry)) + end do + ! y_field <- alpha * A(i_block_row, j_block_col) * x_field + y_field + call a%block_storage%mats(i_block_row, j_block_col)%a%csmv( & + & alpha, x_field, done, y_field, info, trans_op) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='block csmv') + return + end if + ! scatter the row-field output sub-vector back into y + do i_entry = 1, n_owned_row_field + y(a%field_map(i_block_row)%global_local_pos(i_entry)) = y_field(i_entry) + end do + deallocate(y_field) + end if + end do + deallocate(x_field) + end do + + end subroutine psb_d_nest_base_csmv + + ! Selective (regime 2) application of a SINGLE block: + ! y_field = alpha * A(i_block_row, j_block_col) * x_field + beta * y_field + ! x_field is the column-field local vector (owned + ghosts) ALREADY halo-exchanged + ! by the caller; y_field is the row-field owned local vector. The caller chooses + ! the exchange regime (the union halo, or just this block's halo), so this + ! routine is purely local. It is FORMAT-AGNOSTIC: it dispatches to the block's + ! own polymorphic csmv, so the block may be CSR, COO, ... independently of the + ! other blocks. (The full-operator matvec, regime 1, is psb_d_nest_base_csmv.) + subroutine psb_d_nest_apply_block(nest_op, i_block_row, j_block_col, alpha, x_field, beta, y_field, info) + type(psb_d_nest_base_mat), intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: i_block_row, j_block_col + real(psb_dpk_), intent(in) :: alpha, beta, x_field(:) + real(psb_dpk_), intent(inout) :: y_field(:) + integer(psb_ipk_), intent(out) :: info + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_apply_block' + + if (.not. associated(nest_op%block_storage)) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='nested operator not set up') + return + end if + if (.not. nest_op%block_storage%has_block(i_block_row, j_block_col)) then + ! absent block contributes zero: y_field <- beta * y_field + if (beta == dzero) then + y_field(:) = dzero + else if (beta /= done) then + y_field(:) = beta * y_field(:) + end if + return + end if + + ! polymorphic dispatch: the block applies its own format (CSR/COO/...) + call nest_op%block_storage%mats(i_block_row, j_block_col)%a%csmv( & + & alpha, x_field, beta, y_field, info) + if (info /= psb_success_) & + & call psb_errpush(psb_err_from_subroutine_, name, a_err='block csmv') + + end subroutine psb_d_nest_apply_block + + ! ==================================================================== + ! Field-split interface (for the block preconditioner). + ! Exposes the field structure so a fieldsplit/Schur preconditioner can: + ! - know how many fields there are and their owned sizes; + ! - get a block as a standard psb_dspmat_type (sub-preconditioner on A, + ! Schur-complement matvecs with B / B^T); + ! - get a field descriptor (run a field-level Krylov / halo exchange); + ! - restrict the global vector to a field sub-vector and prolong it back. + ! ==================================================================== + + function psb_d_nest_get_n_fields(nest_op) result(n_fields) + type(psb_d_nest_base_mat), intent(in) :: nest_op + integer(psb_ipk_) :: n_fields + n_fields = nest_op%n_fields + end function psb_d_nest_get_n_fields + + function psb_d_nest_get_field_owned(nest_op, field) result(n_owned) + type(psb_d_nest_base_mat), intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: field + integer(psb_ipk_) :: n_owned + n_owned = 0 + if (allocated(nest_op%field_map) .and. field >= 1 .and. field <= nest_op%n_fields) & + & n_owned = nest_op%field_map(field)%n_owned + end function psb_d_nest_get_field_owned + + ! Pointer to block (i,j) as a standard psb_dspmat_type (null if absent). + function psb_d_nest_get_block(nest_op, i_block_row, j_block_col) result(block_ptr) + type(psb_d_nest_base_mat), target, intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: i_block_row, j_block_col + type(psb_dspmat_type), pointer :: block_ptr + block_ptr => null() + if (associated(nest_op%block_storage)) then + if (nest_op%block_storage%has_block(i_block_row, j_block_col)) & + & block_ptr => nest_op%block_storage%mats(i_block_row, j_block_col) + end if + end function psb_d_nest_get_block + + ! Pointer to field k's descriptor (null if not set up). + function psb_d_nest_get_field_desc(nest_op, field) result(desc_ptr) + type(psb_d_nest_base_mat), target, intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: field + type(psb_desc_type), pointer :: desc_ptr + desc_ptr => null() + if (associated(nest_op%grid_desc) .and. field >= 1 .and. field <= nest_op%n_fields) & + & desc_ptr => nest_op%grid_desc%descs(1, field) + end function psb_d_nest_get_field_desc + + ! Restrict: extract field k's OWNED sub-vector from the global local vector. + subroutine psb_d_nest_restrict_field(nest_op, field, x_global, x_field, info) + type(psb_d_nest_base_mat), intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: field + real(psb_dpk_), intent(in) :: x_global(:) + real(psb_dpk_), intent(out) :: x_field(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i_entry, n_owned + info = psb_success_ + if (field < 1 .or. field > nest_op%n_fields) then + info = psb_err_invalid_input_; return + end if + n_owned = nest_op%field_map(field)%n_owned + do i_entry = 1, n_owned + x_field(i_entry) = x_global(nest_op%field_map(field)%global_local_pos(i_entry)) + end do + end subroutine psb_d_nest_restrict_field + + ! Prolong: insert field k's OWNED sub-vector into the global local vector. + subroutine psb_d_nest_prolong_field(nest_op, field, x_field, x_global, info) + type(psb_d_nest_base_mat), intent(in) :: nest_op + integer(psb_ipk_), intent(in) :: field + real(psb_dpk_), intent(in) :: x_field(:) + real(psb_dpk_), intent(inout) :: x_global(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i_entry, n_owned + info = psb_success_ + if (field < 1 .or. field > nest_op%n_fields) then + info = psb_err_invalid_input_; return + end if + n_owned = nest_op%field_map(field)%n_owned + do i_entry = 1, n_owned + x_global(nest_op%field_map(field)%global_local_pos(i_entry)) = x_field(i_entry) + end do + end subroutine psb_d_nest_prolong_field + +end module psb_d_nest_base_mat_mod diff --git a/base/modules/serial/psb_d_nest_mat_mod.f90 b/base/modules/serial/psb_d_nest_mat_mod.f90 index bfcb0c966..dc68513c5 100644 --- a/base/modules/serial/psb_d_nest_mat_mod.f90 +++ b/base/modules/serial/psb_d_nest_mat_mod.f90 @@ -30,13 +30,15 @@ ! ! ! module: psb_d_nest_mat_mod +! Author: Simone Staccone (Stack-1) ! ! Defines psb_d_nest_sparse_mat: a block-structured distributed sparse ! matrix for double precision real arithmetic. ! ! The matrix is stored as a 2-D array of psb_dspmat_type sub-matrices. -! A companion logical array blk_present(i,j) flags which blocks are -! non-null (absent blocks contribute zero to any product). +! Block presence is determined directly from the sub-matrix storage: a block +! (i,j) is present iff mats(i,j)%a is allocated (absent blocks contribute zero +! to any product). There is no separate presence flag array. ! ! Descriptor convention (current nested design) ! --------------------------------------------- @@ -46,7 +48,7 @@ ! descs(i,j) together with mats(i,j). ! ! A block may be structurally absent (NULL/zero): this is represented by -! blk_present(i,j)=.false. and mats(i,j) left unbuilt. In that case the +! mats(i,j) left unbuilt (mats(i,j)%a not allocated). In that case the ! block contributes zero and is skipped by nested kernels. ! ! Descriptor storage is distinct from matrix presence: descriptors are @@ -54,7 +56,7 @@ ! matrix blocks may be present only on a subset. ! ! Reference examples in test/pdegen: -! * psb_d_pde_nest.full.F90 (A(2,2) left NULL, blk_present(2,2)=.false.) +! * psb_d_pde_nest.full.F90 (A(2,2) left NULL, mats(2,2)%a not allocated) ! * psb_d_nest_tools.F90 and psb_d_pde_nest_full_tools.F90 ! (2-D desc_nest%descs(i,j) used in nested allocation/assembly). ! @@ -66,7 +68,6 @@ module psb_d_nest_mat_mod integer(psb_ipk_) :: nrblocks = 0 integer(psb_ipk_) :: ncblocks = 0 type(psb_dspmat_type), allocatable :: mats(:,:) - logical, allocatable :: blk_present(:,:) contains procedure :: get_nrblocks => psb_d_nest_mat_get_nrb procedure :: get_ncblocks => psb_d_nest_mat_get_ncb @@ -91,29 +92,32 @@ contains end function psb_d_nest_mat_get_ncb ! has_block: return .true. if block (i,j) is non-null - function psb_d_nest_mat_has_block(a, i, j) result(hp) + function psb_d_nest_mat_has_block(a, i_block_row, j_block_col) result(has) class(psb_d_nest_sparse_mat), intent(in) :: a - integer(psb_ipk_), intent(in) :: i, j - logical :: hp + integer(psb_ipk_), intent(in) :: i_block_row, j_block_col + logical :: has - hp = .false. - if (i < 1 .or. i > a%nrblocks) return - if (j < 1 .or. j > a%ncblocks) return - if (.not. allocated(a%blk_present)) return - hp = a%blk_present(i, j) + has = .false. + if (i_block_row < 1 .or. i_block_row > a%nrblocks) return + if (j_block_col < 1 .or. j_block_col > a%ncblocks) return + if (.not. allocated(a%mats)) return + ! P3: presence is determined solely by whether the sub-matrix has been + ! built (its polymorphic storage %a is allocated). No parallel flag array. + has = allocated(a%mats(i_block_row, j_block_col)%a) end function psb_d_nest_mat_has_block ! sizeof: total storage across all allocated sub-matrices - function psb_d_nest_mat_sizeof(a) result(s) + function psb_d_nest_mat_sizeof(a) result(total_bytes) class(psb_d_nest_sparse_mat), intent(in) :: a - integer(psb_epk_) :: s - integer(psb_ipk_) :: i, j + integer(psb_epk_) :: total_bytes + integer(psb_ipk_) :: i_block_row, j_block_col - s = 0_psb_epk_ + total_bytes = 0_psb_epk_ if (allocated(a%mats)) then - do j = 1, a%ncblocks - do i = 1, a%nrblocks - if (a%blk_present(i, j)) s = s + a%mats(i, j)%sizeof() + do j_block_col = 1, a%ncblocks + do i_block_row = 1, a%nrblocks + if (allocated(a%mats(i_block_row, j_block_col)%a)) & + & total_bytes = total_bytes + a%mats(i_block_row, j_block_col)%sizeof() end do end do end if @@ -124,23 +128,19 @@ contains class(psb_d_nest_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, j, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, local_info info = 0 if (allocated(a%mats)) then - do j = 1, a%ncblocks - do i = 1, a%nrblocks - if (a%blk_present(i, j)) then - call a%mats(i, j)%free() + do j_block_col = 1, a%ncblocks + do i_block_row = 1, a%nrblocks + if (allocated(a%mats(i_block_row, j_block_col)%a)) then + call a%mats(i_block_row, j_block_col)%free() end if end do end do - deallocate(a%mats, stat=linfo) - if (linfo /= 0 .and. info == 0) info = linfo - end if - if (allocated(a%blk_present)) then - deallocate(a%blk_present, stat=linfo) - if (linfo /= 0 .and. info == 0) info = linfo + deallocate(a%mats, stat=local_info) + if (local_info /= 0 .and. info == 0) info = local_info end if a%nrblocks = 0 a%ncblocks = 0 diff --git a/base/modules/serial/psb_d_nest_vect_mod.f90 b/base/modules/serial/psb_d_nest_vect_mod.f90 deleted file mode 100644 index c90693607..000000000 --- a/base/modules/serial/psb_d_nest_vect_mod.f90 +++ /dev/null @@ -1,109 +0,0 @@ -! -! 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. -! -! -! module: psb_d_nest_vect_mod -! -! Defines psb_d_nest_vect_type: a block-structured distributed dense -! vector for double precision real arithmetic. Each sub-vector is a -! standard psb_d_vect_type assembled under its own descriptor. -! -! Parallel BLAS operations (nrm2, dot, axpby) are exposed as module -! subroutines/functions in psb_d_nest_psblas_mod so that they can -! exploit a single global reduction per call. -! -module psb_d_nest_vect_mod - use psb_d_vect_mod - use psb_desc_mod - implicit none - - type :: psb_d_nest_vect_type - integer(psb_ipk_) :: nblocks = 0 - type(psb_d_vect_type), allocatable :: vects(:) - contains - procedure :: get_nblocks => psb_d_nest_vect_get_nblocks - procedure :: zero => psb_d_nest_vect_zero - procedure :: sizeof => psb_d_nest_vect_sizeof - procedure :: free => psb_d_nest_vect_free - end type psb_d_nest_vect_type - -contains - - ! get_nblocks - function psb_d_nest_vect_get_nblocks(x) result(nb) - class(psb_d_nest_vect_type), intent(in) :: x - integer(psb_ipk_) :: nb - nb = x%nblocks - end function psb_d_nest_vect_get_nblocks - - ! zero: set all sub-vectors to zero (local, no halo zeroing needed) - subroutine psb_d_nest_vect_zero(x) - class(psb_d_nest_vect_type), intent(inout) :: x - integer(psb_ipk_) :: i - if (allocated(x%vects)) then - do i = 1, x%nblocks - call x%vects(i)%zero() - end do - end if - end subroutine psb_d_nest_vect_zero - - ! sizeof: total bytes across all sub-vectors - function psb_d_nest_vect_sizeof(x) result(s) - class(psb_d_nest_vect_type), intent(in) :: x - integer(psb_epk_) :: s - integer(psb_ipk_) :: i - s = 0_psb_epk_ - if (allocated(x%vects)) then - do i = 1, x%nblocks - s = s + x%vects(i)%sizeof() - end do - end if - end function psb_d_nest_vect_sizeof - - ! free: release all sub-vectors - subroutine psb_d_nest_vect_free(x, info) - class(psb_d_nest_vect_type), intent(inout) :: x - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i, linfo - - info = 0 - if (allocated(x%vects)) then - do i = 1, x%nblocks - call x%vects(i)%free(linfo) - if (linfo /= 0 .and. info == 0) info = linfo - end do - deallocate(x%vects, stat=linfo) - if (linfo /= 0 .and. info == 0) info = linfo - end if - x%nblocks = 0 - end subroutine psb_d_nest_vect_free - -end module psb_d_nest_vect_mod diff --git a/base/modules/tools/psb_cd_nest_tools_mod.F90 b/base/modules/tools/psb_cd_nest_tools_mod.F90 index 78454ec32..4fe4487f8 100644 --- a/base/modules/tools/psb_cd_nest_tools_mod.F90 +++ b/base/modules/tools/psb_cd_nest_tools_mod.F90 @@ -29,6 +29,9 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! +! Module: psb_cd_nest_tools_mod +! Author: Simone Staccone (Stack-1) +! ! Nested-specific assembly wrappers for PSBLAS3 — descriptor routines ! @@ -39,12 +42,14 @@ module psb_cd_nest_tools_mod use psb_error_mod, only : psb_errpush use psb_cd_tools_mod, only : psb_cdall, psb_cdasb, psb_cdins, psb_cdcpy, psb_cdprt use psb_desc_nest_mod, only : psb_desc_nest_type + use psb_desc_mod, only : psb_desc_type implicit none private public :: psb_cdall_nest, psb_cdins_nest, psb_cdins_nest_rc, & - psb_cdasb_nest, psb_cdfree_nest, psb_cdcpy_nest, psb_cdprt_nest + psb_cdasb_nest, psb_cdfree_nest, psb_cdcpy_nest, psb_cdprt_nest, & + psb_cd_nest_compose ! Column-only form: (blk_j, nz, ja, desc_nest, info [,mask, lidx]) ! Row+column form: (blk_i, blk_j, nz, ia, ja, desc_nest, info) @@ -86,29 +91,29 @@ contains integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nrblocks, ncblocks, nl - integer(psb_ipk_) :: i, j, nr, nc, nl_ + integer(psb_ipk_) :: i_block_row, j_block_col, n_block_rows, n_block_cols, local_rows character(len=20) :: name info = psb_success_ name = 'psb_cdall_nest' ! Set default dimensions - nr = 2 - nc = 2 - if (present(nrblocks)) nr = nrblocks - if (present(ncblocks)) nc = ncblocks + n_block_rows = 2 + n_block_cols = 2 + if (present(nrblocks)) n_block_rows = nrblocks + if (present(ncblocks)) n_block_cols = ncblocks if (.not. present(nl)) then info = psb_err_no_optional_arg_ call psb_errpush(info, name, a_err='nl (local row count)') return end if - nl_ = nl + local_rows = nl ! Allocate nested descriptor structure - desc_nest%nrblocks = nr - desc_nest%ncblocks = nc - allocate(desc_nest%descs(nr, nc), stat=info) + desc_nest%nrblocks = n_block_rows + desc_nest%ncblocks = n_block_cols + allocate(desc_nest%descs(n_block_rows, n_block_cols), stat=info) if (info /= 0) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) @@ -120,9 +125,9 @@ contains ! psb_cdasb_nest assembles both the original and the clone the shared ! base_desc is rebuilt twice, corrupting the global-to-local mapping of ! every block in that row. Independent allocations avoid this entirely. - do i = 1, nr - do j = 1, nc - call psb_cdall(ctxt, desc_nest%descs(i, j), info, nl=nl_) + do i_block_row = 1, n_block_rows + do j_block_col = 1, n_block_cols + call psb_cdall(ctxt, desc_nest%descs(i_block_row, j_block_col), info, nl=local_rows) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_, name) return @@ -135,9 +140,9 @@ contains #if defined(PSB_IPK4) && defined(PSB_LPK8) ! psb_cdins_nest_rc_sub: row+col form, ipk_ nz — only when ipk_ /= lpk_ - subroutine psb_cdins_nest_rc_sub(blk_i, blk_j, nz, ia, ja, desc_nest, info) - integer(psb_ipk_), intent(in) :: blk_i, blk_j, nz - integer(psb_lpk_), intent(in) :: ia(:), ja(:) + subroutine psb_cdins_nest_rc_sub(block_row, block_col, n_entries, entry_rows, entry_cols, desc_nest, info) + integer(psb_ipk_), intent(in) :: block_row, block_col, n_entries + integer(psb_lpk_), intent(in) :: entry_rows(:), entry_cols(:) type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info @@ -146,16 +151,16 @@ contains info = psb_success_ name = 'psb_cdins_nest' - if (nz == 0) return + if (n_entries == 0) return - if (blk_i < 1 .or. blk_i > desc_nest%nrblocks .or. & - blk_j < 1 .or. blk_j > desc_nest%ncblocks) then + if (block_row < 1 .or. block_row > desc_nest%nrblocks .or. & + block_col < 1 .or. block_col > desc_nest%ncblocks) then info = psb_err_invalid_input_ call psb_errpush(info, name, a_err='invalid block indices') return end if - call psb_cdins(nz, ia, ja, desc_nest%descs(blk_i, blk_j), info) + call psb_cdins(n_entries, entry_rows, entry_cols, desc_nest%descs(block_row, block_col), info) if (info /= psb_success_) & call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins') @@ -163,45 +168,45 @@ contains ! psb_cdins_nest_c: col-only form, ipk_ nz — only when ipk_ /= lpk_ - subroutine psb_cdins_nest_c(blk_j, nz, ja, desc_nest, info, mask, lidx) - integer(psb_ipk_), intent(in) :: blk_j, nz - integer(psb_lpk_), intent(in) :: ja(:) + subroutine psb_cdins_nest_c(block_col, n_entries, entry_cols, desc_nest, info, mask, lidx) + integer(psb_ipk_), intent(in) :: block_col, n_entries + integer(psb_lpk_), intent(in) :: entry_cols(:) type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info logical, intent(in), optional, target :: mask(:) integer(psb_ipk_), intent(in), optional :: lidx(:) - integer(psb_ipk_) :: i, linfo + integer(psb_ipk_) :: i_block_row, local_info character(len=20) :: name info = psb_success_ name = 'psb_cdins_nest' - if (nz == 0) return + if (n_entries == 0) return - if (blk_j < 1 .or. blk_j > desc_nest%ncblocks) then + if (block_col < 1 .or. block_col > desc_nest%ncblocks) then info = psb_err_invalid_input_ call psb_errpush(info, name, a_err='invalid block column index') return end if - do i = 1, desc_nest%nrblocks - linfo = psb_success_ + do i_block_row = 1, desc_nest%nrblocks + local_info = psb_success_ if (present(mask)) then if (present(lidx)) then - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, mask=mask, lidx=lidx) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, mask=mask, lidx=lidx) else - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, mask=mask) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, mask=mask) end if else if (present(lidx)) then - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, lidx=lidx) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, lidx=lidx) else - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info) end if end if - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo + if (local_info /= psb_success_ .and. info == psb_success_) then + info = local_info call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins') end if end do @@ -215,9 +220,9 @@ contains ! When entries in block (blk_i, blk_j) reference columns owned by other ! processes, use the col-only form afterwards to broadcast those column ! indices across all row-blocks in block-col blk_j. - subroutine psb_lcdins_nest_rc(blk_i, blk_j, nz, ia, ja, desc_nest, info) - integer(psb_ipk_), intent(in) :: blk_i, blk_j - integer(psb_lpk_), intent(in) :: nz, ia(:), ja(:) + subroutine psb_lcdins_nest_rc(block_row, block_col, n_entries, entry_rows, entry_cols, desc_nest, info) + integer(psb_ipk_), intent(in) :: block_row, block_col + integer(psb_lpk_), intent(in) :: n_entries, entry_rows(:), entry_cols(:) type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info @@ -226,16 +231,16 @@ contains info = psb_success_ name = 'psb_cdins_nest' - if (nz == 0) return + if (n_entries == 0) return - if (blk_i < 1 .or. blk_i > desc_nest%nrblocks .or. & - blk_j < 1 .or. blk_j > desc_nest%ncblocks) then + if (block_row < 1 .or. block_row > desc_nest%nrblocks .or. & + block_col < 1 .or. block_col > desc_nest%ncblocks) then info = psb_err_invalid_input_ call psb_errpush(info, name, a_err='invalid block indices') return end if - call psb_cdins(nz, ia, ja, desc_nest%descs(blk_i, blk_j), info) + call psb_cdins(n_entries, entry_rows, entry_cols, desc_nest%descs(block_row, block_col), info) if (info /= psb_success_) & call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins') @@ -247,45 +252,45 @@ contains ! Registers nz global column indices ja into the descriptor for ! block column blk_j across all row-blocks (descs(i, blk_j) for ! i = 1..nrblocks). mask and lidx are forwarded to psb_cdins. - subroutine psb_lcdins_nest_c(blk_j, nz, ja, desc_nest, info, mask, lidx) - integer(psb_ipk_), intent(in) :: blk_j - integer(psb_lpk_), intent(in) :: nz, ja(:) + subroutine psb_lcdins_nest_c(block_col, n_entries, entry_cols, desc_nest, info, mask, lidx) + integer(psb_ipk_), intent(in) :: block_col + integer(psb_lpk_), intent(in) :: n_entries, entry_cols(:) type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info logical, intent(in), optional, target :: mask(:) integer(psb_ipk_), intent(in), optional :: lidx(:) - integer(psb_ipk_) :: i, linfo + integer(psb_ipk_) :: i_block_row, local_info character(len=20) :: name info = psb_success_ name = 'psb_cdins_nest' - if (nz == 0) return + if (n_entries == 0) return - if (blk_j < 1 .or. blk_j > desc_nest%ncblocks) then + if (block_col < 1 .or. block_col > desc_nest%ncblocks) then info = psb_err_invalid_input_ call psb_errpush(info, name, a_err='invalid block column index') return end if - do i = 1, desc_nest%nrblocks - linfo = psb_success_ + do i_block_row = 1, desc_nest%nrblocks + local_info = psb_success_ if (present(mask)) then if (present(lidx)) then - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, mask=mask, lidx=lidx) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, mask=mask, lidx=lidx) else - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, mask=mask) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, mask=mask) end if else if (present(lidx)) then - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo, lidx=lidx) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info, lidx=lidx) else - call psb_cdins(nz, ja, desc_nest%descs(i, blk_j), linfo) + call psb_cdins(n_entries, entry_cols, desc_nest%descs(i_block_row, block_col), local_info) end if end if - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo + if (local_info /= psb_success_ .and. info == psb_success_) then + info = local_info call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins') end if end do @@ -307,15 +312,15 @@ contains type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: i_block_row, j_block_col character(len=20) :: name info = psb_success_ name = 'psb_cdasb_nest' - do i = 1, desc_nest%nrblocks - do j = 1, desc_nest%ncblocks - call psb_cdasb(desc_nest%descs(i, j), info) + do i_block_row = 1, desc_nest%nrblocks + do j_block_col = 1, desc_nest%ncblocks + call psb_cdasb(desc_nest%descs(i_block_row, j_block_col), info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdasb') return @@ -368,7 +373,7 @@ contains type(psb_desc_nest_type), intent(out) :: desc_out integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: i_block_row, j_block_col character(len=20) :: name info = psb_success_ @@ -383,9 +388,9 @@ contains return end if - do i = 1, desc_in%nrblocks - do j = 1, desc_in%ncblocks - call psb_cdcpy(desc_in%descs(i, j), desc_out%descs(i, j), info) + do i_block_row = 1, desc_in%nrblocks + do j_block_col = 1, desc_in%ncblocks + call psb_cdcpy(desc_in%descs(i_block_row, j_block_col), desc_out%descs(i_block_row, j_block_col), info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdcpy') return @@ -415,37 +420,37 @@ contains logical, intent(in), optional :: glob, short integer(psb_ipk_), intent(in), optional :: verbosity - integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: i_block_row, j_block_col - do i = 1, desc_nest%nrblocks - do j = 1, desc_nest%ncblocks - write(iout, '(a,i0,a,i0,a)') 'Block (', i, ',', j, '):' + do i_block_row = 1, desc_nest%nrblocks + do j_block_col = 1, desc_nest%ncblocks + write(iout, '(a,i0,a,i0,a)') 'Block (', i_block_row, ',', j_block_col, '):' if (present(glob)) then if (present(short)) then if (present(verbosity)) then - call psb_cdprt(iout, desc_nest%descs(i,j), glob=glob, short=short, verbosity=verbosity) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), glob=glob, short=short, verbosity=verbosity) else - call psb_cdprt(iout, desc_nest%descs(i,j), glob=glob, short=short) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), glob=glob, short=short) end if else if (present(verbosity)) then - call psb_cdprt(iout, desc_nest%descs(i,j), glob=glob, verbosity=verbosity) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), glob=glob, verbosity=verbosity) else - call psb_cdprt(iout, desc_nest%descs(i,j), glob=glob) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), glob=glob) end if end if else if (present(short)) then if (present(verbosity)) then - call psb_cdprt(iout, desc_nest%descs(i,j), short=short, verbosity=verbosity) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), short=short, verbosity=verbosity) else - call psb_cdprt(iout, desc_nest%descs(i,j), short=short) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), short=short) end if else if (present(verbosity)) then - call psb_cdprt(iout, desc_nest%descs(i,j), verbosity=verbosity) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col), verbosity=verbosity) else - call psb_cdprt(iout, desc_nest%descs(i,j)) + call psb_cdprt(iout, desc_nest%descs(i_block_row,j_block_col)) end if end if end if @@ -454,4 +459,126 @@ contains end subroutine psb_cdprt_nest + ! psb_cd_nest_compose (P1 / step 6a) + ! + ! Compose the per-field block descriptors into a SINGLE global psb_desc_type + ! describing the whole nested operator. The global index space is the + ! concatenation of the field spaces: + ! + ! global index = offset_k + (field-k global index), offset_k = sum_{m 0) then + allocate(halo_global(halo_count), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + halo_count = 0 + do i_field = 1, n_fields + n_owned = desc_grid%descs(1,i_field)%get_local_rows() + n_local = desc_grid%descs(1,i_field)%get_local_cols() + do i_loc = n_owned + 1, n_local + call desc_grid%descs(1,i_field)%l2g(i_loc, global_idx, info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='l2g halo'); return + end if + halo_count = halo_count + 1 + halo_global(halo_count) = field_offset(i_field) + global_idx + end do + end do + call psb_cdins(halo_count, halo_global, desc_global, info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins'); return + end if + end if + + ! 5. assemble: build the global halo communication schedule (union halo) + call psb_cdasb(desc_global, info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdasb'); return + end if + + end subroutine psb_cd_nest_compose + end module psb_cd_nest_tools_mod diff --git a/base/modules/tools/psb_d_nest_builder_mod.F90 b/base/modules/tools/psb_d_nest_builder_mod.F90 new file mode 100644 index 000000000..0fb2851ea --- /dev/null +++ b/base/modules/tools/psb_d_nest_builder_mod.F90 @@ -0,0 +1,368 @@ +! +! 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 prior 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_nest_builder_mod.F90 +! +! Module: psb_d_nest_builder_mod +! Author: Simone Staccone (Stack-1) +! +! User-friendly frontend to build a nested (MATNEST) operator without manually +! managing per-field descriptors, the union halo, composition and setup. +! +! All the boilerplate (identical for every nested operator) is hidden behind a +! single type, psb_d_nest_matrix, with the usual PSBLAS init/ins/asb pattern: +! +! type(psb_d_nest_matrix) :: nested_matrix +! call nested_matrix%init(ctxt, [n1, n2], info) ! 2 fields of global size n1, n2 +! call nested_matrix%ins(1,1, n, rows, cols, vals, info) ! values of block (1,1) = A +! call nested_matrix%ins(1,2, n, rows, cols, vals, info) ! values of block (1,2) = B^T +! call nested_matrix%ins(2,1, n, rows, cols, vals, info) ! values of block (2,1) = B +! ... ! (absent blocks = not inserted) +! call nested_matrix%asb(info) ! assemble: builds a_glob, desc_glob +! +! ! from here on nested_matrix%a_glob and nested_matrix%desc_glob are an +! ! ordinary distributed matrix/descriptor: +! call psb_geall(x, nested_matrix%desc_glob, info) +! call psb_krylov('CG', nested_matrix%a_glob, prec, b, x, eps, nested_matrix%desc_glob, info, ...) +! +! Indices: in ins(block_row, block_col, ...) the rows live in the index space of +! field block_row, the columns in the index space of field block_col (GLOBAL +! field indices, 1..field_size). Each process inserts only the rows it owns +! (PSBLAS convention). Off-diagonal blocks may be rectangular. +! +! NOTE: after asb the object holds consistent internal pointers (a_glob%a points +! to block_storage / grid_desc): do not copy/move the object after assembly. +! +module psb_d_nest_builder_mod + use psb_const_mod + use psb_error_mod, only : psb_errpush + use psb_penv_mod, only : psb_ctxt_type, psb_info + use psb_desc_mod, only : psb_desc_type + use psb_d_mat_mod, only : psb_dspmat_type + use psb_cd_tools_mod, only : psb_cdall, psb_cdins, psb_cdasb + use psb_desc_nest_mod, only : psb_desc_nest_type + use psb_d_nest_mat_mod, only : psb_d_nest_sparse_mat + use psb_d_nest_base_mat_mod, only : psb_d_nest_base_mat, psb_d_nest_base_setup + use psb_cd_nest_tools_mod, only : psb_cd_nest_compose + use psb_d_nest_tools_mod, only : psb_d_nest_rect_block + implicit none + + ! growing triplet buffer for a single block + type :: psb_d_nest_block_buffer + integer(psb_ipk_) :: n_entries = 0 + integer(psb_lpk_), allocatable :: entry_rows(:), entry_cols(:) + real(psb_dpk_), allocatable :: entry_vals(:) + end type psb_d_nest_block_buffer + + type :: psb_d_nest_matrix + type(psb_ctxt_type) :: context + integer(psb_ipk_) :: n_fields = 0 + logical :: assembled = .false. + ! construction state + type(psb_desc_type), allocatable :: field_desc(:) ! one descriptor per field + type(psb_d_nest_block_buffer), allocatable :: block_buffer(:,:) ! triplets per block (i,j) + ! products (owned; the pointers in a_glob%a point in here) + type(psb_d_nest_sparse_mat) :: block_storage + type(psb_desc_nest_type) :: grid_desc + type(psb_dspmat_type) :: a_glob ! the matrix to hand to Krylov + type(psb_desc_type) :: desc_glob ! the global descriptor + contains + procedure, pass(op) :: init => psb_d_nest_op_init + procedure, pass(op) :: ins => psb_d_nest_op_ins + procedure, pass(op) :: asb => psb_d_nest_op_asb + procedure, pass(op) :: free => psb_d_nest_op_free + end type psb_d_nest_matrix + + private + public :: psb_d_nest_matrix + +contains + + ! init: create one descriptor per field (block distribution from the global sizes) + subroutine psb_d_nest_op_init(op, context, field_sizes, info) + class(psb_d_nest_matrix), intent(inout) :: op + type(psb_ctxt_type), intent(in) :: context + integer(psb_lpk_), intent(in) :: field_sizes(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: my_rank, num_procs, n_fields, i_field, field_local_rows + integer(psb_lpk_) :: field_global_size + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_op_init' + + call psb_info(context, my_rank, num_procs) + n_fields = size(field_sizes) + op%context = context + op%n_fields = n_fields + op%assembled = .false. + + allocate(op%field_desc(n_fields), op%block_buffer(n_fields,n_fields), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + + do i_field = 1, n_fields + field_global_size = field_sizes(i_field) + ! block distribution: field_global_size rows over num_procs processes (total size invariant) + field_local_rows = int(field_global_size / int(num_procs, psb_lpk_), psb_ipk_) + if (int(my_rank, psb_lpk_) < mod(field_global_size, int(num_procs, psb_lpk_))) & + & field_local_rows = field_local_rows + 1 + call psb_cdall(context, op%field_desc(i_field), info, nl=field_local_rows) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdall'); return + end if + end do + end subroutine psb_d_nest_op_init + + ! ins: accumulate the triplets into block (block_row,block_col) and register the + ! columns (field block_col index space) into that descriptor's union halo + subroutine psb_d_nest_op_ins(op, block_row, block_col, n_entries, entry_rows, entry_cols, entry_vals, info) + class(psb_d_nest_matrix), intent(inout) :: op + integer(psb_ipk_), intent(in) :: block_row, block_col, n_entries + integer(psb_lpk_), intent(in) :: entry_rows(:), entry_cols(:) + real(psb_dpk_), intent(in) :: entry_vals(:) + integer(psb_ipk_), intent(out) :: info + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_op_ins' + + if (op%assembled) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='operator already assembled'); return + end if + if (block_row < 1 .or. block_row > op%n_fields .or. & + & block_col < 1 .or. block_col > op%n_fields) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='block index out of range'); return + end if + if (n_entries <= 0) return + + call block_buffer_append(op%block_buffer(block_row,block_col), n_entries, & + & entry_rows, entry_cols, entry_vals, info) + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + + ! the columns of block (block_row,block_col) live in field block_col -> + ! register their indices into that descriptor's union halo + ! (this also applies when block_col == block_row) + call psb_cdins(n_entries, entry_cols(1:n_entries), op%field_desc(block_col), info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdins'); return + end if + end subroutine psb_d_nest_op_ins + + ! asb: assemble the descriptors, build the blocks, compose the global + ! descriptor, set up the operator and wrap it into a_glob + subroutine psb_d_nest_op_asb(op, info) + class(psb_d_nest_matrix), intent(inout), target :: op + integer(psb_ipk_), intent(out) :: info + + type(psb_d_nest_base_mat) :: nest_operator + integer(psb_ipk_) :: n_fields, i_field, j_field + character(len=24) :: name + + info = psb_success_ + name = 'psb_d_nest_op_asb' + n_fields = op%n_fields + + ! 1) assemble the per-field descriptors (with the union halo accumulated in ins) + do i_field = 1, n_fields + call psb_cdasb(op%field_desc(i_field), info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdasb'); return + end if + end do + + ! 2) build the local blocks (generally rectangular) from the triplets + op%block_storage%nrblocks = n_fields + op%block_storage%ncblocks = n_fields + allocate(op%block_storage%mats(n_fields,n_fields), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + do j_field = 1, n_fields + do i_field = 1, n_fields + if (op%block_buffer(i_field,j_field)%n_entries > 0) then + call psb_d_nest_rect_block(op%block_storage%mats(i_field,j_field), & + & op%block_buffer(i_field,j_field)%n_entries, & + & op%block_buffer(i_field,j_field)%entry_rows, & + & op%block_buffer(i_field,j_field)%entry_cols, & + & op%block_buffer(i_field,j_field)%entry_vals, & + & op%field_desc(i_field), op%field_desc(j_field), info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='rect_block'); return + end if + end if + end do + end do + + ! 3) descriptor grid: descs(i,j) = descriptor of field j + op%grid_desc%nrblocks = n_fields + op%grid_desc%ncblocks = n_fields + allocate(op%grid_desc%descs(n_fields,n_fields), stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + do j_field = 1, n_fields + do i_field = 1, n_fields + call op%field_desc(j_field)%clone(op%grid_desc%descs(i_field,j_field), info) + end do + end do + + ! 4) composed global descriptor + operator setup + call psb_cd_nest_compose(op%grid_desc, op%desc_glob, info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='cd_nest_compose'); return + end if + call psb_d_nest_base_setup(nest_operator, op%block_storage, op%grid_desc, op%desc_glob, info) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='nest_base_setup'); return + end if + + ! 5) wrap into the standard matrix object (the pointers keep pointing at op%*) + allocate(op%a_glob%a, source=nest_operator, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return + end if + call op%a_glob%set_nrows(op%desc_glob%get_local_rows()) + call op%a_glob%set_ncols(op%desc_glob%get_local_cols()) + call op%a_glob%set_asb() + + ! 6) the triplet buffers are no longer needed + do j_field = 1, n_fields + do i_field = 1, n_fields + call block_buffer_free(op%block_buffer(i_field,j_field)) + end do + end do + op%assembled = .true. + end subroutine psb_d_nest_op_asb + + ! free: release everything + subroutine psb_d_nest_op_free(op, info) + class(psb_d_nest_matrix), intent(inout) :: op + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i_field, j_field, local_info + + info = psb_success_ + if (allocated(op%block_buffer)) then + do j_field = 1, size(op%block_buffer,2) + do i_field = 1, size(op%block_buffer,1) + call block_buffer_free(op%block_buffer(i_field,j_field)) + end do + end do + deallocate(op%block_buffer, stat=local_info) + end if + if (op%assembled) then + call op%a_glob%free() + call op%desc_glob%free(local_info) + call op%grid_desc%free(local_info) + end if + if (allocated(op%field_desc)) then + do i_field = 1, size(op%field_desc) + call op%field_desc(i_field)%free(local_info) + end do + deallocate(op%field_desc, stat=local_info) + end if + op%n_fields = 0 + op%assembled = .false. + end subroutine psb_d_nest_op_free + + !----------------------------------------------------------------- + ! private helpers: growing triplet buffer + !----------------------------------------------------------------- + subroutine block_buffer_append(buffer, n_entries, entry_rows, entry_cols, entry_vals, info) + type(psb_d_nest_block_buffer), intent(inout) :: buffer + integer(psb_ipk_), intent(in) :: n_entries + integer(psb_lpk_), intent(in) :: entry_rows(:), entry_cols(:) + real(psb_dpk_), intent(in) :: entry_vals(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: required_size + + info = psb_success_ + required_size = buffer%n_entries + n_entries + call ensure_capacity_lpk(buffer%entry_rows, required_size, info); if (info /= 0) return + call ensure_capacity_lpk(buffer%entry_cols, required_size, info); if (info /= 0) return + call ensure_capacity_dpk(buffer%entry_vals, required_size, info); if (info /= 0) return + buffer%entry_rows(buffer%n_entries+1:required_size) = entry_rows(1:n_entries) + buffer%entry_cols(buffer%n_entries+1:required_size) = entry_cols(1:n_entries) + buffer%entry_vals(buffer%n_entries+1:required_size) = entry_vals(1:n_entries) + buffer%n_entries = required_size + end subroutine block_buffer_append + + subroutine ensure_capacity_lpk(array, required_size, info) + integer(psb_lpk_), allocatable, intent(inout) :: array(:) + integer(psb_ipk_), intent(in) :: required_size + integer(psb_ipk_), intent(out) :: info + integer(psb_lpk_), allocatable :: grown(:) + integer(psb_ipk_) :: capacity + + info = 0 + if (.not. allocated(array)) then + allocate(array(max(required_size,16)), stat=info); return + end if + capacity = size(array) + if (required_size <= capacity) return + allocate(grown(max(2*capacity, required_size)), stat=info); if (info /= 0) return + grown(1:capacity) = array(1:capacity) + call move_alloc(grown, array) + end subroutine ensure_capacity_lpk + + subroutine ensure_capacity_dpk(array, required_size, info) + real(psb_dpk_), allocatable, intent(inout) :: array(:) + integer(psb_ipk_), intent(in) :: required_size + integer(psb_ipk_), intent(out) :: info + real(psb_dpk_), allocatable :: grown(:) + integer(psb_ipk_) :: capacity + + info = 0 + if (.not. allocated(array)) then + allocate(array(max(required_size,16)), stat=info); return + end if + capacity = size(array) + if (required_size <= capacity) return + allocate(grown(max(2*capacity, required_size)), stat=info); if (info /= 0) return + grown(1:capacity) = array(1:capacity) + call move_alloc(grown, array) + end subroutine ensure_capacity_dpk + + subroutine block_buffer_free(buffer) + type(psb_d_nest_block_buffer), intent(inout) :: buffer + if (allocated(buffer%entry_rows)) deallocate(buffer%entry_rows) + if (allocated(buffer%entry_cols)) deallocate(buffer%entry_cols) + if (allocated(buffer%entry_vals)) deallocate(buffer%entry_vals) + buffer%n_entries = 0 + end subroutine block_buffer_free + +end module psb_d_nest_builder_mod diff --git a/base/modules/tools/psb_d_nest_tools_mod.F90 b/base/modules/tools/psb_d_nest_tools_mod.F90 index 687a46c2d..f1e8f2836 100644 --- a/base/modules/tools/psb_d_nest_tools_mod.F90 +++ b/base/modules/tools/psb_d_nest_tools_mod.F90 @@ -29,6 +29,9 @@ ! POSSIBILITY OF SUCH DAMAGE. ! ! +! Module: psb_d_nest_tools_mod +! Author: Simone Staccone (Stack-1) +! ! Nested-specific assembly wrappers for PSBLAS3 — double precision matrix and vector routines ! @@ -41,13 +44,15 @@ module psb_d_nest_tools_mod psb_geall, psb_geins, psb_geasb, psb_gefree use psb_desc_nest_mod, only : psb_desc_nest_type use psb_d_nest_mat_mod, only : psb_d_nest_sparse_mat - use psb_d_nest_vect_mod, only : psb_d_nest_vect_type + use psb_d_mat_mod, only : psb_dspmat_type + use psb_d_base_mat_mod, only : psb_d_coo_sparse_mat + use psb_desc_mod, only : psb_desc_type implicit none private public :: psb_spall_nest, psb_spins_nest, psb_spasb_nest, psb_spfree_nest, psb_sprn_nest, & - psb_geall_nest, psb_geins_nest, psb_geasb_nest, psb_gefree_nest + psb_d_nest_rect_block contains @@ -62,7 +67,7 @@ contains integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nnz - integer(psb_ipk_) :: i, j, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, local_info character(len=20) :: name info = psb_success_ @@ -80,30 +85,19 @@ contains end if end if - if (.not. allocated(a_nest%blk_present)) then - allocate(a_nest%blk_present(a_nest%nrblocks, a_nest%ncblocks), stat=info) - if (info /= 0) then - info = psb_err_alloc_dealloc_ - call psb_errpush(info, name) - return - end if - a_nest%blk_present = .false. - end if - - do i = 1, a_nest%nrblocks - do j = 1, a_nest%ncblocks - linfo = psb_success_ + do i_block_row = 1, a_nest%nrblocks + do j_block_col = 1, a_nest%ncblocks + local_info = psb_success_ if (present(nnz)) then - call psb_spall(a_nest%mats(i, j), desc_nest%descs(i, j), linfo, nnz=nnz) + call psb_spall(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), local_info, nnz=nnz) else - call psb_spall(a_nest%mats(i, j), desc_nest%descs(i, j), linfo) + call psb_spall(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), local_info) end if - if (linfo /= psb_success_) then - info = linfo + if (local_info /= psb_success_) then + info = local_info call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spall') return end if - a_nest%blk_present(i, j) = .true. end do end do @@ -114,24 +108,24 @@ contains ! The block is lazy-allocated on first insertion if psb_spall_nest ! was not called first. - subroutine psb_spins_nest(blk_i, blk_j, nz, ia, ja, val, a_nest, desc_nest, info) - integer(psb_ipk_), intent(in) :: blk_i, blk_j, nz - integer(psb_lpk_), intent(in) :: ia(:), ja(:) - real(psb_dpk_), intent(in) :: val(:) + subroutine psb_spins_nest(block_row, block_col, n_entries, entry_rows, entry_cols, entry_vals, a_nest, desc_nest, info) + integer(psb_ipk_), intent(in) :: block_row, block_col, n_entries + integer(psb_lpk_), intent(in) :: entry_rows(:), entry_cols(:) + real(psb_dpk_), intent(in) :: entry_vals(:) type(psb_d_nest_sparse_mat), intent(inout) :: a_nest type(psb_desc_nest_type), intent(inout) :: desc_nest integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: nnz_est + integer(psb_ipk_) :: nnz_estimate character(len=20) :: name info = psb_success_ name = 'psb_spins_nest' - if (nz == 0) return + if (n_entries == 0) return - if (blk_i < 1 .or. blk_i > a_nest%nrblocks .or. & - blk_j < 1 .or. blk_j > a_nest%ncblocks) then + if (block_row < 1 .or. block_row > a_nest%nrblocks .or. & + block_col < 1 .or. block_col > a_nest%ncblocks) then info = psb_err_invalid_input_ call psb_errpush(info, name, a_err='invalid block indices') return @@ -144,29 +138,21 @@ contains call psb_errpush(info, name) return end if - allocate(a_nest%blk_present(a_nest%nrblocks, a_nest%ncblocks), stat=info) - if (info /= 0) then - info = psb_err_alloc_dealloc_ - call psb_errpush(info, name) - return - end if - a_nest%blk_present = .false. end if - if (.not. a_nest%blk_present(blk_i, blk_j)) then - ! Estimate nnz: use nz + 50% buffer for future insertions - nnz_est = max(nz, 10) + nz / 2 - call psb_spall(a_nest%mats(blk_i, blk_j), & - desc_nest%descs(blk_i, blk_j), info, nnz=nnz_est) + if (.not. allocated(a_nest%mats(block_row, block_col)%a)) then + ! Estimate nnz: use n_entries + 50% buffer for future insertions + nnz_estimate = max(n_entries, 10) + n_entries / 2 + call psb_spall(a_nest%mats(block_row, block_col), & + desc_nest%descs(block_row, block_col), info, nnz=nnz_estimate) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spall') return end if - a_nest%blk_present(blk_i, blk_j) = .true. end if - call psb_spins(nz, ia, ja, val, a_nest%mats(blk_i, blk_j), & - desc_nest%descs(blk_i, blk_j), info) + call psb_spins(n_entries, entry_rows, entry_cols, entry_vals, a_nest%mats(block_row, block_col), & + desc_nest%descs(block_row, block_col), info) if (info /= psb_success_) & call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spins') @@ -181,32 +167,33 @@ contains integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: dupl - integer(psb_ipk_) :: i, j, dupl_, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, dupl_mode, local_info character(len=20) :: name - info = psb_success_ - name = 'psb_spasb_nest' - dupl_ = psb_dupl_add_ - if (present(dupl)) dupl_ = dupl - - do i = 1, a_nest%nrblocks - do j = 1, a_nest%ncblocks - if (a_nest%blk_present(i, j)) then - linfo = psb_success_ - if (dupl_ == psb_dupl_add_) then - call psb_spasb(a_nest%mats(i, j), desc_nest%descs(i, j), linfo, & - dupl=psb_dupl_add_) - else if (dupl_ == psb_dupl_ovwrt_) then - call psb_spasb(a_nest%mats(i, j), desc_nest%descs(i, j), linfo, & - dupl=psb_dupl_ovwrt_) - else if (dupl_ == psb_dupl_err_) then - call psb_spasb(a_nest%mats(i, j), desc_nest%descs(i, j), linfo, & - dupl=psb_dupl_err_) + info = psb_success_ + name = 'psb_spasb_nest' + dupl_mode = psb_dupl_add_ + if (present(dupl)) dupl_mode = dupl + + do i_block_row = 1, a_nest%nrblocks + do j_block_col = 1, a_nest%ncblocks + if (allocated(a_nest%mats(i_block_row, j_block_col)%a)) then + local_info = psb_success_ + if (dupl_mode == psb_dupl_add_) then + call psb_spasb(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), & + local_info, dupl=psb_dupl_add_) + else if (dupl_mode == psb_dupl_ovwrt_) then + call psb_spasb(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), & + local_info, dupl=psb_dupl_ovwrt_) + else if (dupl_mode == psb_dupl_err_) then + call psb_spasb(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), & + local_info, dupl=psb_dupl_err_) else - call psb_spasb(a_nest%mats(i, j), desc_nest%descs(i, j), linfo) + call psb_spasb(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), & + local_info) end if - if (linfo /= psb_success_) then - info = linfo + if (local_info /= psb_success_) then + info = local_info call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spasb') return end if @@ -217,44 +204,34 @@ contains end subroutine psb_spasb_nest ! Calls psb_spfree on every present block, then deallocates the - ! mats and blk_present arrays and resets nrblocks/ncblocks to 0. + ! mats array and resets nrblocks/ncblocks to 0. subroutine psb_spfree_nest(a_nest, desc_nest, info) type(psb_d_nest_sparse_mat), intent(inout) :: a_nest type(psb_desc_nest_type), intent(in) :: desc_nest integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, j, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, local_info character(len=20) :: name info = psb_success_ name = 'psb_spfree_nest' if (allocated(a_nest%mats)) then - do i = 1, a_nest%nrblocks - do j = 1, a_nest%ncblocks - if (allocated(a_nest%blk_present)) then - if (a_nest%blk_present(i, j)) then - linfo = psb_success_ - call psb_spfree(a_nest%mats(i, j), desc_nest%descs(i, j), linfo) - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo - call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spfree') - end if + do i_block_row = 1, a_nest%nrblocks + do j_block_col = 1, a_nest%ncblocks + if (allocated(a_nest%mats(i_block_row, j_block_col)%a)) then + local_info = psb_success_ + call psb_spfree(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), local_info) + if (local_info /= psb_success_ .and. info == psb_success_) then + info = local_info + call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_spfree') end if end if end do end do - deallocate(a_nest%mats, stat=linfo) - if (linfo /= 0 .and. info == psb_success_) then - info = psb_err_alloc_dealloc_ - call psb_errpush(info, name) - end if - end if - - if (allocated(a_nest%blk_present)) then - deallocate(a_nest%blk_present, stat=linfo) - if (linfo /= 0 .and. info == psb_success_) then + deallocate(a_nest%mats, stat=local_info) + if (local_info /= 0 .and. info == psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) end if @@ -274,25 +251,25 @@ contains integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: clear - integer(psb_ipk_) :: i, j, linfo + integer(psb_ipk_) :: i_block_row, j_block_col, local_info character(len=20) :: name info = psb_success_ name = 'psb_sprn_nest' - if (.not. allocated(a_nest%mats) .or. .not. allocated(a_nest%blk_present)) return + if (.not. allocated(a_nest%mats)) return - do i = 1, a_nest%nrblocks - do j = 1, a_nest%ncblocks - if (a_nest%blk_present(i, j)) then - linfo = psb_success_ + do i_block_row = 1, a_nest%nrblocks + do j_block_col = 1, a_nest%ncblocks + if (allocated(a_nest%mats(i_block_row, j_block_col)%a)) then + local_info = psb_success_ if (present(clear)) then - call psb_sprn(a_nest%mats(i, j), desc_nest%descs(i, j), linfo, clear=clear) + call psb_sprn(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), local_info, clear=clear) else - call psb_sprn(a_nest%mats(i, j), desc_nest%descs(i, j), linfo) + call psb_sprn(a_nest%mats(i_block_row, j_block_col), desc_nest%descs(i_block_row, j_block_col), local_info) end if - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo + if (local_info /= psb_success_ .and. info == psb_success_) then + info = local_info call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_sprn') end if end if @@ -302,132 +279,78 @@ contains end subroutine psb_sprn_nest - ! Allocates one sub-vector per block-row, using descs(i, 1) as - ! the row descriptor for block i. Must be called before psb_geins_nest. - - subroutine psb_geall_nest(x_nest, desc_nest, info) - type(psb_d_nest_vect_type), intent(out) :: x_nest - type(psb_desc_nest_type), intent(in) :: desc_nest - integer(psb_ipk_), intent(out) :: info - integer(psb_ipk_) :: i, linfo - character(len=20) :: name + ! psb_d_nest_rect_block + ! + ! Build a local GENERAL (possibly rectangular) block A(i,j) of a nested + ! operator, with rows in field i and columns in field j (field i /= field j, + ! |field i| /= |field j| allowed). Rows are localized against the field-i + ! (row) descriptor, columns against the field-j (column) descriptor — which + ! must already carry the union halo of column j (cdall + cdins(all column-j + ! blocks' columns) + cdasb). The result is a CSR block of shape + ! (field-i owned rows) x (field-j local cols incl. halo) + ! consumable directly by the nested csmv (psb_d_nest_base_mat). + ! + ! A single-descriptor psb_spall/psb_spasb cannot express row-field /= col-field + ! (it would force rows and columns into the same index space), hence the + ! explicit COO build with separate row/column localization. + ! + ! Arguments (this process's local contribution): + ! blk (out) the assembled block (CSR) + ! nz number of local entries + ! ia_glob(:) GLOBAL field-i row indices (owned by this process) + ! ja_glob(:) GLOBAL field-j column indices + ! val(:) values + ! desc_row field-i descriptor (rows) + ! desc_col field-j descriptor (columns, with union halo) + ! + subroutine psb_d_nest_rect_block(blk, nz, ia_glob, ja_glob, val, desc_row, desc_col, info) + type(psb_dspmat_type), intent(out) :: blk + integer(psb_ipk_), intent(in) :: nz + integer(psb_lpk_), intent(in) :: ia_glob(:), ja_glob(:) + real(psb_dpk_), intent(in) :: val(:) + type(psb_desc_type), intent(in) :: desc_row, desc_col + integer(psb_ipk_), intent(out) :: info + + type(psb_d_coo_sparse_mat) :: coo_block + integer(psb_ipk_) :: k_entry, n_loc_rows, n_loc_cols, loc_row, loc_col + character(len=24) :: name info = psb_success_ - name = 'psb_geall_nest' + name = 'psb_d_nest_rect_block' - x_nest%nblocks = desc_nest%nrblocks - allocate(x_nest%vects(x_nest%nblocks), stat=info) - if (info /= 0) then - info = psb_err_alloc_dealloc_ - call psb_errpush(info, name) - return - end if + n_loc_rows = desc_row%get_local_rows() ! owned rows of field i + n_loc_cols = desc_col%get_local_cols() ! field-j local cols (owned + halo) - do i = 1, x_nest%nblocks - linfo = psb_success_ - call psb_geall(x_nest%vects(i), desc_nest%descs(i, 1), linfo) - if (linfo /= psb_success_) then - info = linfo - call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_geall') + call coo_block%allocate(n_loc_rows, n_loc_cols, nz) + do k_entry = 1, nz + call desc_row%g2l(ia_glob(k_entry), loc_row, info) + if (info /= 0 .or. loc_row < 1 .or. loc_row > n_loc_rows) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='row not owned / not localizable') return end if - end do - - end subroutine psb_geall_nest - - ! Inserts m entries into block blk_i of the nested vector. - - subroutine psb_geins_nest(blk_i, m, irw, val, x_nest, desc_nest, info, local) - integer(psb_ipk_), intent(in) :: blk_i, m - integer(psb_lpk_), intent(in) :: irw(:) - real(psb_dpk_), intent(in) :: val(:) - type(psb_d_nest_vect_type), intent(inout) :: x_nest - type(psb_desc_nest_type), intent(in) :: desc_nest - integer(psb_ipk_), intent(out) :: info - logical, intent(in), optional :: local - - character(len=20) :: name - - info = psb_success_ - name = 'psb_geins_nest' - - if (m == 0) return - - if (blk_i < 1 .or. blk_i > x_nest%nblocks) then - info = psb_err_invalid_input_ - call psb_errpush(info, name, a_err='invalid block index') - return - end if - - if (present(local)) then - call psb_geins(m, irw, val, x_nest%vects(blk_i), desc_nest%descs(blk_i, 1), info, & - local=local) - else - call psb_geins(m, irw, val, x_nest%vects(blk_i), desc_nest%descs(blk_i, 1), info) - end if - if (info /= psb_success_) & - call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_geins') - - end subroutine psb_geins_nest - - ! Calls psb_geasb on every sub-vector. - ! Must be called after psb_cdasb_nest and all psb_geins_nest calls. - - subroutine psb_geasb_nest(x_nest, desc_nest, info) - type(psb_d_nest_vect_type), intent(inout) :: x_nest - type(psb_desc_nest_type), intent(in) :: desc_nest - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i, linfo - character(len=20) :: name - - info = psb_success_ - name = 'psb_geasb_nest' - - do i = 1, x_nest%nblocks - linfo = psb_success_ - call psb_geasb(x_nest%vects(i), desc_nest%descs(i, 1), linfo) - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo - call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_geasb') + call desc_col%g2l(ja_glob(k_entry), loc_col, info) + if (info /= 0 .or. loc_col < 1 .or. loc_col > n_loc_cols) then + info = psb_err_invalid_input_ + call psb_errpush(info, name, a_err='column not in field-j descriptor (missing from union halo)') + return end if + coo_block%ia(k_entry) = loc_row + coo_block%ja(k_entry) = loc_col + coo_block%val(k_entry) = val(k_entry) end do - - end subroutine psb_geasb_nest - - ! Calls psb_gefree on every sub-vector, then deallocates the - ! vects array and resets nblocks to 0. - - subroutine psb_gefree_nest(x_nest, desc_nest, info) - type(psb_d_nest_vect_type), intent(inout) :: x_nest - type(psb_desc_nest_type), intent(in) :: desc_nest - integer(psb_ipk_), intent(out) :: info - - integer(psb_ipk_) :: i, linfo - character(len=20) :: name - - info = psb_success_ - name = 'psb_gefree_nest' - - if (allocated(x_nest%vects)) then - do i = 1, x_nest%nblocks - linfo = psb_success_ - call psb_gefree(x_nest%vects(i), desc_nest%descs(i, 1), linfo) - if (linfo /= psb_success_ .and. info == psb_success_) then - info = linfo - call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_gefree') - end if - end do - deallocate(x_nest%vects, stat=linfo) - if (linfo /= 0 .and. info == psb_success_) then - info = psb_err_alloc_dealloc_ - call psb_errpush(info, name) - end if + call coo_block%set_nzeros(nz) + call coo_block%set_dupl(psb_dupl_add_) + call coo_block%fix(info) + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='coo fix'); return end if - - x_nest%nblocks = 0 - - end subroutine psb_gefree_nest + call blk%mv_from(coo_block) + call blk%cscnv(info, type='CSR') + if (info /= 0) then + call psb_errpush(psb_err_from_subroutine_, name, a_err='cscnv'); return + end if + end subroutine psb_d_nest_rect_block end module psb_d_nest_tools_mod diff --git a/test/nested/CMakeLists.txt b/test/nested/CMakeLists.txt new file mode 100644 index 000000000..95c1bed75 --- /dev/null +++ b/test/nested/CMakeLists.txt @@ -0,0 +1,48 @@ +cmake_minimum_required(VERSION 3.10) +project(nested Fortran) + +# Check for the installation path for psblas +if(NOT DEFINED PSBLAS_INSTALL_DIR) + message(FATAL_ERROR "Please specify the path to the psblas installation directory using -DPSBLAS_INSTALL_DIR=") +endif() + +# Set the include and library directories based on the provided path +set(INSTALLDIR "${PSBLAS_INSTALL_DIR}") +set(INCDIR "${INSTALLDIR}/include") +set(MODDIR "${INSTALLDIR}/modules") +set(LIBDIR "${INSTALLDIR}/lib") + +# Find the psblas package +find_package(psblas REQUIRED PATHS ${INSTALLDIR}) + +# Include directories for the Fortran compiler +include_directories(${INCDIR} ${MODDIR}) + +# Define executable directory +set(EXEDIR "${CMAKE_CURRENT_SOURCE_DIR}/runs") +file(MAKE_DIRECTORY ${EXEDIR}) + +# Nested (block-structured / MATNEST) tests +set(SOURCES_D_NEST_GLOB_TEST psb_d_nest_glob_test.F90) +set(SOURCES_D_NEST_RECT_TEST psb_d_nest_rect_test.F90) +set(SOURCES_D_NEST_CG_TEST psb_d_nest_cg_test.F90) +set(SOURCES_D_NEST_BUILDER_TEST psb_d_nest_builder_test.F90) + +add_executable(psb_d_nest_glob_test ${SOURCES_D_NEST_GLOB_TEST}) +target_link_libraries(psb_d_nest_glob_test psblas::util psblas::linsolve psblas::prec psblas::base) + +add_executable(psb_d_nest_rect_test ${SOURCES_D_NEST_RECT_TEST}) +target_link_libraries(psb_d_nest_rect_test psblas::util psblas::linsolve psblas::prec psblas::base) + +add_executable(psb_d_nest_cg_test ${SOURCES_D_NEST_CG_TEST}) +target_link_libraries(psb_d_nest_cg_test psblas::util psblas::linsolve psblas::prec psblas::base) + +add_executable(psb_d_nest_builder_test ${SOURCES_D_NEST_BUILDER_TEST}) +target_link_libraries(psb_d_nest_builder_test psblas::util psblas::linsolve psblas::prec psblas::base) + +# Set output directory for executables +foreach(target psb_d_nest_glob_test psb_d_nest_rect_test psb_d_nest_cg_test psb_d_nest_builder_test) + set_target_properties(${target} PROPERTIES + RUNTIME_OUTPUT_DIRECTORY ${EXEDIR} + ) +endforeach() diff --git a/test/nested/Makefile b/test/nested/Makefile new file mode 100644 index 000000000..b59501ecd --- /dev/null +++ b/test/nested/Makefile @@ -0,0 +1,46 @@ +INSTALLDIR=../.. +INCDIR=$(INSTALLDIR)/include +MODDIR=$(INSTALLDIR)/modules/ +include $(INCDIR)/Make.inc.psblas +# +# Libraries used +LIBDIR=$(INSTALLDIR)/lib +PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_linsolve -lpsb_prec -lpsb_base +LDLIBS= $(PSBLDLIBS) +# +# Compilers and such +# +CCOPT= -g +FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG). + + +EXEDIR=./runs + +all: runsd psb_d_nest_glob_test psb_d_nest_rect_test psb_d_nest_cg_test psb_d_nest_builder_test + +runsd: + (if test ! -d runs ; then mkdir runs; fi) + +psb_d_nest_glob_test: psb_d_nest_glob_test.o + $(FLINK) psb_d_nest_glob_test.o -o psb_d_nest_glob_test $(PSBLAS_LIB) $(LDLIBS) + /bin/mv psb_d_nest_glob_test $(EXEDIR) + +psb_d_nest_rect_test: psb_d_nest_rect_test.o + $(FLINK) psb_d_nest_rect_test.o -o psb_d_nest_rect_test $(PSBLAS_LIB) $(LDLIBS) + /bin/mv psb_d_nest_rect_test $(EXEDIR) + +psb_d_nest_cg_test: psb_d_nest_cg_test.o + $(FLINK) psb_d_nest_cg_test.o -o psb_d_nest_cg_test $(PSBLAS_LIB) $(LDLIBS) + /bin/mv psb_d_nest_cg_test $(EXEDIR) + +psb_d_nest_builder_test: psb_d_nest_builder_test.o + $(FLINK) psb_d_nest_builder_test.o -o psb_d_nest_builder_test $(PSBLAS_LIB) $(LDLIBS) + /bin/mv psb_d_nest_builder_test $(EXEDIR) + +clean: + /bin/rm -f psb_d_nest_glob_test.o psb_d_nest_rect_test.o psb_d_nest_cg_test.o psb_d_nest_builder_test.o *$(.mod) \ + $(EXEDIR)/psb_d_nest_glob_test $(EXEDIR)/psb_d_nest_rect_test $(EXEDIR)/psb_d_nest_cg_test $(EXEDIR)/psb_d_nest_builder_test +verycleanlib: + (cd ../..; make veryclean) +lib: + (cd ../../; make library) diff --git a/test/nested/README.md b/test/nested/README.md new file mode 100644 index 000000000..a52d2cc5c --- /dev/null +++ b/test/nested/README.md @@ -0,0 +1,133 @@ +# Nested (block-structured / MATNEST) matrices in PSBLAS + +Author: Simone Staccone (Stack-1) + +This directory contains the tests for the **nested matrix** support added to PSBLAS: a block-structured distributed operator + +``` + [ A11 A12 ... ] + M = [ A21 A22 ... ] + [ ... ... ] +``` + +whose blocks are kept as separate sparse matrices (one per field) but which presents itself to Krylov solvers and preconditioners as a **single ordinary distributed matrix**. It is the PSBLAS analogue of PETSc's `MATNEST`. + +The motivating case is the **saddle-point** system + +``` + M = [ A B^T ] + [ B 0 ] +``` + +(symmetric indefinite, with the (2,2) block absent), but the implementation supports any square multi-field block operator with possibly **rectangular** +sub-blocks. + + +## 1. Concepts + +* **Field** — a contiguous index space (e.g. velocity `V` and pressure `Q` in a saddle-point problem). Each field has its own `psb_desc_type` distribution. +* **Block (i,j)** — the sub-matrix coupling field `i` (rows) with field `j` (columns). It may be rectangular (`|field i| /= |field j|`) and may be absent. +* **Global operator** — the blocks are concatenated into a single **square** operator `M` of size `sum(field_sizes)`, distributed over one **composed global descriptor** with a **union halo** (one halo exchange per matrix-vector product, covering all blocks of a given column field at once). +* **Rectangular blocks** — PSBLAS does not support rectangular *distributed* matrices, but it does support rectangular *local* CSR/COO matrices. The rectangular product therefore happens only in the **local** block `csmv`; the only object carrying a descriptor (and hence communication) is the global operator, which is always square. + +The global operator (`a_glob`) and global descriptor (`desc_glob`) can be passed unchanged to `psb_spmm`, `psb_krylov`, and the standard preconditioners. + + +## 2. Recommended API: `psb_d_nest_matrix` + +The easy way to build a nested matrix is the `psb_d_nest_matrix` type (module `psb_d_nest_builder_mod`, re-exported by the umbrella `psb_d_nest_mod`), which follows the usual PSBLAS `init` / `ins` / `asb` pattern and hides all the descriptor / halo / compose / setup boilerplate: + +```fortran +use psb_d_nest_mod + +type(psb_d_nest_matrix) :: nested_matrix +integer(psb_lpk_) :: n1, n2 + +! 1) declare the field structure: 2 fields of global size n1, n2 +call nested_matrix%init(ctxt, [n1, n2], info) + +! 2) insert the block values, owned rows only (PSBLAS convention). +! ins(block_row, block_col, n_entries, entry_rows, entry_cols, entry_vals, info) +! rows are GLOBAL indices in field block_row, columns in field block_col. +call nested_matrix%ins(1, 1, nz_A, iaA, jaA, valA, info) ! A = block (1,1) +call nested_matrix%ins(1, 2, nz_Bt, iaBt, jaBt, valBt, info) ! B^T = block (1,2) +call nested_matrix%ins(2, 1, nz_B, iaB, jaB, valB, info) ! B = block (2,1) +! (the (2,2) block is simply not inserted) + +! 3) assemble: builds nested_matrix%a_glob and nested_matrix%desc_glob +call nested_matrix%asb(info) + +! 4) from here on it is an ordinary distributed matrix/descriptor +call psb_geall(x, nested_matrix%desc_glob, info) +... +call psb_krylov('CG', nested_matrix%a_glob, prec, b, x, eps, & + & nested_matrix%desc_glob, info, itmax=..., iter=..., err=...) + +! 5) release +call nested_matrix%free(info) +``` + +Notes: + +* To know which rows it owns in a field, a process can query the per-field descriptor exposed as `nested_matrix%field_desc(i)` (e.g. `nested_matrix%field_desc(1)%get_local_rows()` and `%l2g(...)`), exactly as it would with a plain `psb_cdall` descriptor. +* Off-diagonal blocks may be rectangular: the cross-field column indices are registered into the union halo automatically by `ins`. +* The CG solver requires an SPD operator; a genuine saddle-point operator is indefinite and needs MINRES/GMRES (plus, eventually, a block preconditioner). +* **Do not copy/move** a `psb_d_nest_matrix` after `asb`: the wrapped operator holds internal pointers into the object. + + +## 3. Low-level path (advanced) + +`psb_d_nest_matrix` is built on three lower-level pieces, available directly for advanced use (see `psb_d_nest_cg_test.F90` for an end-to-end example): + +* `psb_cd_nest_compose(grid_desc, desc_glob, info)` — compose the per-field descriptors into the single global descriptor with the union halo. +* `psb_d_nest_base_setup(nest_op, block_storage, grid_desc, desc_glob, info)` — set up the `psb_d_nest_base_mat` operator (implements the local `csmv`). +* `psb_d_nest_rect_block(blk, nz, ia, ja, val, desc_row, desc_col, info)` — build a single (possibly rectangular) local block from global triplets, with rows localized against `desc_row` and columns against `desc_col`. + +A field-split interface (`psb_d_nest_get_block`, `psb_d_nest_get_field_desc`, +`psb_d_nest_restrict_field`, `psb_d_nest_prolong_field`, +`psb_d_nest_apply_block`) is exposed on `psb_d_nest_base_mat` as the hook for a future block (field-split / Schur) preconditioner. + + +## 4. Tests + +| Test | What it checks | +|------------------------------|----------------| +| `psb_d_nest_glob_test` | Square 2×2 operator built with `psb_d_nest_matrix`; the nested `psb_spmm` is compared bit-for-bit against the same matrix assembled monolithically in CSR. | +| `psb_d_nest_rect_test` | Same, with fields of different size (`|V| = 2|Q|`) and genuinely **rectangular** off-diagonal blocks. | +| `psb_d_nest_cg_test` | Standard PSBLAS **CG** on an SPD, ill-conditioned operator (1D Laplacian reordered red-black), built on the **low-level path**; the solution is recovered to machine precision over hundreds of matvecs. | +| `psb_d_nest_builder_test` | Same CG solve as above but built through the `psb_d_nest_matrix` utility (high-level path). | + +All tests run both serially and in parallel, and the result is invariant with respect to the number of MPI processes. + +### Build and run + +The PSBLAS library must be built/installed first (from the repository root): + +```sh +make # or the CMake build +``` + +Then, from this directory: + +```sh +make # builds the executables into ./runs +./runs/psb_d_nest_glob_test # serial +mpirun -np 4 ./runs/psb_d_nest_rect_test +mpirun -np 4 ./runs/psb_d_nest_cg_test +mpirun -np 4 ./runs/psb_d_nest_builder_test +``` + +Each test prints a single `[PASS]` / `[FAIL]` line (printed by rank 0). + + +## 5. Source files + +Library (under `base/modules/`): + +* `desc/psb_desc_nest_mod.f90` — `psb_desc_nest_type` (grid of per-field descriptors) +* `serial/psb_d_nest_mat_mod.f90` — `psb_d_nest_sparse_mat` (block storage) +* `serial/psb_d_nest_base_mat_mod.F90`— `psb_d_nest_base_mat` (the MATNEST operator + `csmv`) +* `tools/psb_cd_nest_tools_mod.F90` — descriptor tools (`psb_cd_nest_compose`, ...) +* `tools/psb_d_nest_tools_mod.F90` — block tools (`psb_d_nest_rect_block`, ...) +* `tools/psb_d_nest_builder_mod.F90` — `psb_d_nest_matrix` frontend (init/ins/asb) +* `psb_d_nest_mod.f90` — umbrella module (`use psb_d_nest_mod`) diff --git a/test/nested/psb_d_nest_builder_test.F90 b/test/nested/psb_d_nest_builder_test.F90 new file mode 100644 index 000000000..9a03aceff --- /dev/null +++ b/test/nested/psb_d_nest_builder_test.F90 @@ -0,0 +1,207 @@ +! +! 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 prior 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_nest_builder_test.F90 +! +! Program: psb_d_nest_builder_test +! Author: Simone Staccone (Stack-1) +! +! Same operator as the low-level CG test (1D Laplacian reordered red-black, SPD +! and ill-conditioned) but built with the psb_d_nest_matrix utility: the user +! declares nested_matrix, gives the field sizes, inserts the block values and +! calls asb. All the setup (per-field descriptors, union halo, compose, setup, +! wrap) is handled by the utility. Solved with CG and checked against the +! exact solution. +! +! M = [ 2I C ] C(r,r) = -1 , C(r,r-1) = -1 (the Laplacian edges) +! [ C^T 2I ] +! +! Run: ./psb_d_nest_builder_test ; mpirun -np 4 ./psb_d_nest_builder_test +! +program psb_d_nest_builder_test + use psb_base_mod + use psb_prec_mod + use psb_linsolve_mod + use psb_d_nest_mod ! umbrella: includes psb_d_nest_matrix (builder) + implicit none + + type(psb_ctxt_type) :: context + integer(psb_ipk_) :: my_rank, num_procs, info, i_local_row, entry_idx + integer(psb_ipk_) :: field1_local_rows, field2_local_rows + integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size + + type(psb_d_nest_matrix) :: nested_matrix ! the only object needed + type(psb_dprec_type) :: preconditioner + type(psb_d_vect_type) :: x_solution, rhs, x_exact + real(psb_dpk_) :: insert_value(1) + + integer(psb_lpk_), allocatable :: entry_rows(:), entry_cols(:) + real(psb_dpk_), allocatable :: entry_vals(:) + + real(psb_dpk_) :: stop_tol, final_residual, norm_x_exact, solution_error + integer(psb_ipk_) :: max_iter, n_iter, stop_criterion + real(psb_dpk_), parameter :: solution_tol = 1.0e-6_psb_dpk_ + + call psb_init(context) + call psb_info(context, my_rank, num_procs) + + field_size = 512 ! global size of each field (N = 2*field_size) + stop_tol = 1.0e-9_psb_dpk_ + max_iter = 4000 + stop_criterion = 2 + + !--------------------------------------------------------------- + ! 1) create the nested operator: 2 fields of global size field_size + !--------------------------------------------------------------- + call nested_matrix%init(context, [field_size, field_size], info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL init info=', info; goto 9999 + end if + + ! rows owned by this process in each field + field1_local_rows = nested_matrix%field_desc(1)%get_local_rows() + field2_local_rows = nested_matrix%field_desc(2)%get_local_rows() + + !--------------------------------------------------------------- + ! 2) insert the values, one block at a time (owned rows only) + !--------------------------------------------------------------- + ! block (1,1) = 2I + allocate(entry_rows(field1_local_rows), entry_cols(field1_local_rows), entry_vals(field1_local_rows)) + do i_local_row = 1, field1_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, field1_global_row, info) + entry_rows(i_local_row)=field1_global_row; entry_cols(i_local_row)=field1_global_row + entry_vals(i_local_row)=2.0_psb_dpk_ + end do + call nested_matrix%ins(1, 1, field1_local_rows, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! block (2,2) = 2I + allocate(entry_rows(field2_local_rows), entry_cols(field2_local_rows), entry_vals(field2_local_rows)) + do i_local_row = 1, field2_local_rows + call nested_matrix%field_desc(2)%l2g(i_local_row, field2_global_row, info) + entry_rows(i_local_row)=field2_global_row; entry_cols(i_local_row)=field2_global_row + entry_vals(i_local_row)=2.0_psb_dpk_ + end do + call nested_matrix%ins(2, 2, field2_local_rows, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! block (1,2) = C : rows field1, cols field2 ; C(r,r)=-1, C(r,r-1)=-1 + allocate(entry_rows(2*field1_local_rows), entry_cols(2*field1_local_rows), entry_vals(2*field1_local_rows)) + entry_idx = 0 + do i_local_row = 1, field1_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, field1_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field1_global_row + entry_cols(entry_idx) = field1_global_row + entry_vals(entry_idx) = -1.0_psb_dpk_ + if (field1_global_row > 1) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field1_global_row + entry_cols(entry_idx) = field1_global_row - 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! block (2,1) = C^T : rows field2, cols field1 ; C^T(s,s)=-1, C^T(s,s+1)=-1 + allocate(entry_rows(2*field2_local_rows), entry_cols(2*field2_local_rows), entry_vals(2*field2_local_rows)) + entry_idx = 0 + do i_local_row = 1, field2_local_rows + call nested_matrix%field_desc(2)%l2g(i_local_row, field2_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field2_global_row + entry_cols(entry_idx) = field2_global_row + entry_vals(entry_idx) = -1.0_psb_dpk_ + if (field2_global_row < field_size) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field2_global_row + entry_cols(entry_idx) = field2_global_row + 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + !--------------------------------------------------------------- + ! 3) assemble: from here nested_matrix%a_glob / nested_matrix%desc_glob are ready for Krylov + !--------------------------------------------------------------- + call nested_matrix%asb(info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL asb info=', info; goto 9999 + end if + + !--------------------------------------------------------------- + ! 4) consistent RHS x_exact=1, rhs = M*x_exact, then solve with standard CG + !--------------------------------------------------------------- + call psb_geall(x_exact, nested_matrix%desc_glob, info) + do i_local_row = 1, nested_matrix%desc_glob%get_local_rows() + call nested_matrix%desc_glob%l2g(i_local_row, field1_global_row, info) + insert_value(1) = 1.0_psb_dpk_ + call psb_geins(1, [field1_global_row], insert_value, x_exact, nested_matrix%desc_glob, info) + end do + call psb_geasb(x_exact, nested_matrix%desc_glob, info) + + call psb_geall(rhs, nested_matrix%desc_glob, info); call psb_geasb(rhs, nested_matrix%desc_glob, info) + call psb_spmm(done, nested_matrix%a_glob, x_exact, dzero, rhs, nested_matrix%desc_glob, info) + norm_x_exact = psb_genrm2(x_exact, nested_matrix%desc_glob, info) + + call preconditioner%init(context, 'NONE', info) + call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info) + + call psb_geall(x_solution, nested_matrix%desc_glob, info); call psb_geasb(x_solution, nested_matrix%desc_glob, info) + call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, nested_matrix%desc_glob, info, & + & itmax=max_iter, iter=n_iter, err=final_residual, istop=stop_criterion) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL krylov info=', info; goto 9999 + end if + + call psb_geaxpby(-done, x_exact, done, x_solution, nested_matrix%desc_glob, info) + solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info) / norm_x_exact + + if (my_rank == 0) then + write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(global)=', 2*field_size + write(*,'(a,i0)') ' CG iterations = ', n_iter + write(*,'(a,es12.4)') ' CG relative residual = ', final_residual + write(*,'(a,es12.4)') ' ||x - x_exact||/||x_ex|| = ', solution_error + if ((n_iter < max_iter) .and. (solution_error <= solution_tol)) then + write(*,*) '[PASS] nested matrix built with the utility, solved with CG' + else + write(*,*) '[FAIL] tol ', solution_tol + end if + end if + + call nested_matrix%free(info) + +9999 continue + call psb_exit(context) + +end program psb_d_nest_builder_test diff --git a/test/nested/psb_d_nest_cg_test.F90 b/test/nested/psb_d_nest_cg_test.F90 new file mode 100644 index 000000000..d49098fae --- /dev/null +++ b/test/nested/psb_d_nest_cg_test.F90 @@ -0,0 +1,301 @@ +! +! 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 prior 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_nest_cg_test.F90 +! +! Program: psb_d_nest_cg_test +! Author: Simone Staccone (Stack-1) +! +! Solves a linear system with the GLOBAL nested operator using the standard +! PSBLAS CG (psb_krylov('CG', ...)). This test builds the operator on the +! LOW-LEVEL path (per-field descriptors, blocks, compose, setup, wrap) to +! directly validate the machinery the psb_d_nest_matrix utility relies on; the +! same solve through the utility is in psb_d_nest_builder_test. +! +! CG needs a SYMMETRIC POSITIVE DEFINITE operator and, to stress the test +! (hundreds of matvecs), an ILL-CONDITIONED one. We use a real case: the 1D +! Laplacian tridiag(-1, 2, -1) on m = 2*field_size nodes, REORDERED red-black +! (odd nodes -> field 1, even nodes -> field 2). Under this reordering the +! Laplacian becomes exactly +! +! M = [ 2I C ] C(r,r) = -1 , C(r,r-1) = -1 (the Laplacian edges) +! [ C^T 2I ] C^T = exact transpose +! +! (odd nodes are not adjacent to each other -> diagonal blocks = 2I; every -1 +! edge of the Laplacian becomes the coupling C). M is therefore the 1D +! Laplacian up to a permutation: SPD but with lambda_min ~ (pi/m)^2 => cond ~ +! N^2 => CG performs O(N) iterations that GROW with N. +! +! Run: ./psb_d_nest_cg_test ; mpirun -np 4 ./psb_d_nest_cg_test +! +program psb_d_nest_cg_test + use psb_base_mod + use psb_util_mod + use psb_prec_mod + use psb_linsolve_mod + use psb_d_nest_mod + implicit none + + type(psb_ctxt_type) :: context + integer(psb_ipk_) :: my_rank, num_procs, info, i_local_row, entry_idx, field_local_rows + integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size + + ! per-field descriptors + blocks + type(psb_desc_type) :: field1_desc, field2_desc + type(psb_dspmat_type) :: diag_block1, coupling_12, coupling_21, diag_block2 + + ! nested storage + grid descriptor + composed global path + type(psb_d_nest_sparse_mat) :: block_storage + type(psb_desc_nest_type) :: grid_desc + type(psb_desc_type) :: desc_global + type(psb_d_nest_base_mat) :: nest_operator + type(psb_dspmat_type) :: global_operator + + ! preconditioner + vectors + type(psb_dprec_type) :: preconditioner + type(psb_d_vect_type) :: x_solution, rhs, x_exact + real(psb_dpk_) :: insert_value(1) + + ! global triplets for the coupling blocks + integer(psb_lpk_), allocatable :: entry_rows(:), entry_cols(:) + real(psb_dpk_), allocatable :: entry_vals(:) + + ! solver parameters + real(psb_dpk_) :: diag_value, stop_tol, final_residual, norm_x_exact, solution_error + integer(psb_ipk_) :: max_iter, trace_level, n_iter, stop_criterion + real(psb_dpk_), parameter :: solution_tol = 1.0e-6_psb_dpk_ + + call psb_init(context) + call psb_info(context, my_rank, num_procs) + + field_size = 512 ! global rows per field (global N = 2*field_size) + diag_value = 2.0_psb_dpk_ ! Laplacian diagonal (diagonal blocks = diag*I) + stop_tol = 1.0e-9_psb_dpk_ + max_iter = 4000 + trace_level = 0 + stop_criterion = 2 ! stop on the relative residual + + !--------------------------------------------------------------- + ! 1) per-field descriptors: block distribution of field_size global rows + ! over num_procs processes (total size independent of num_procs) + !--------------------------------------------------------------- + field_local_rows = int(field_size / int(num_procs, psb_lpk_), psb_ipk_) + if (int(my_rank, psb_lpk_) < mod(field_size, int(num_procs, psb_lpk_))) & + & field_local_rows = field_local_rows + 1 + call psb_cdall(context, field1_desc, info, nl=field_local_rows) + call psb_cdall(context, field2_desc, info, nl=field_local_rows) + + !--------------------------------------------------------------- + ! 2) diagonal blocks A = B = diag*I (odd/even nodes of the red-black + ! reordered Laplacian are not adjacent to each other) + !--------------------------------------------------------------- + call psb_spall(diag_block1, field1_desc, info, nnz=field1_desc%get_local_rows()) + call psb_spall(diag_block2, field2_desc, info, nnz=field2_desc%get_local_rows()) + + do i_local_row = 1, field1_desc%get_local_rows() + call field1_desc%l2g(i_local_row, field1_global_row, info) + insert_value(1) = diag_value + call psb_spins(1,[field1_global_row],[field1_global_row],insert_value,diag_block1,field1_desc,info) + end do + do i_local_row = 1, field2_desc%get_local_rows() + call field2_desc%l2g(i_local_row, field2_global_row, info) + insert_value(1) = diag_value + call psb_spins(1,[field2_global_row],[field2_global_row],insert_value,diag_block2,field2_desc,info) + end do + + !--------------------------------------------------------------- + ! 3) register, in the union halo, the cross-field columns of the coupling blocks + ! C (row field1, col field2): columns {r, r-1} in field2 -> into field2_desc + ! C^T (row field2, col field1): columns {s, s+1} in field1 -> into field1_desc + !--------------------------------------------------------------- + do i_local_row = 1, field1_desc%get_local_rows() + call field1_desc%l2g(i_local_row, field1_global_row, info) + call psb_cdins(1, [field1_global_row], field2_desc, info) + if (field1_global_row > 1) call psb_cdins(1, [field1_global_row-1_psb_lpk_], field2_desc, info) + end do + do i_local_row = 1, field2_desc%get_local_rows() + call field2_desc%l2g(i_local_row, field2_global_row, info) + call psb_cdins(1, [field2_global_row], field1_desc, info) + if (field2_global_row < field_size) call psb_cdins(1, [field2_global_row+1_psb_lpk_], field1_desc, info) + end do + + call psb_cdasb(field1_desc, info) + call psb_cdasb(field2_desc, info) + call psb_spasb(diag_block1, field1_desc, info, dupl=psb_dupl_add_) + call psb_spasb(diag_block2, field2_desc, info, dupl=psb_dupl_add_) + + !--------------------------------------------------------------- + ! 4) coupling C (1,2): rows field1 (field1_desc), columns field2 (field2_desc) + ! C(r,r) = -1 , C(r,r-1) = -1 (odd node 2r-1 -> even nodes 2r and 2r-2) + !--------------------------------------------------------------- + allocate(entry_rows(2*field1_desc%get_local_rows()), entry_cols(2*field1_desc%get_local_rows()), & + & entry_vals(2*field1_desc%get_local_rows())) + entry_idx = 0 + do i_local_row = 1, field1_desc%get_local_rows() + call field1_desc%l2g(i_local_row, field1_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field1_global_row + entry_cols(entry_idx) = field1_global_row + entry_vals(entry_idx) = -1.0_psb_dpk_ + if (field1_global_row > 1) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field1_global_row + entry_cols(entry_idx) = field1_global_row - 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call psb_d_nest_rect_block(coupling_12, entry_idx, entry_rows, entry_cols, entry_vals, field1_desc, field2_desc, info) + deallocate(entry_rows, entry_cols, entry_vals) + + !--------------------------------------------------------------- + ! 5) coupling C^T (2,1) = exact transpose of C: + ! rows field2 (field2_desc), columns field1 (field1_desc) + ! C^T(s,s) = -1 , C^T(s,s+1) = -1 (even node 2s -> odd nodes 2s-1 and 2s+1) + !--------------------------------------------------------------- + allocate(entry_rows(2*field2_desc%get_local_rows()), entry_cols(2*field2_desc%get_local_rows()), & + & entry_vals(2*field2_desc%get_local_rows())) + entry_idx = 0 + do i_local_row = 1, field2_desc%get_local_rows() + call field2_desc%l2g(i_local_row, field2_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field2_global_row + entry_cols(entry_idx) = field2_global_row + entry_vals(entry_idx) = -1.0_psb_dpk_ + if (field2_global_row < field_size) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = field2_global_row + entry_cols(entry_idx) = field2_global_row + 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call psb_d_nest_rect_block(coupling_21, entry_idx, entry_rows, entry_cols, entry_vals, field2_desc, field1_desc, info) + deallocate(entry_rows, entry_cols, entry_vals) + + !--------------------------------------------------------------- + ! 6) nested grid (all four blocks present) + !--------------------------------------------------------------- + block_storage%nrblocks = 2 + block_storage%ncblocks = 2 + allocate(block_storage%mats(2,2)) + call psb_move_alloc(diag_block1, block_storage%mats(1,1), info) + call psb_move_alloc(coupling_12, block_storage%mats(1,2), info) + call psb_move_alloc(coupling_21, block_storage%mats(2,1), info) + call psb_move_alloc(diag_block2, block_storage%mats(2,2), info) + + grid_desc%nrblocks = 2 + grid_desc%ncblocks = 2 + allocate(grid_desc%descs(2,2)) + call field1_desc%clone(grid_desc%descs(1,1), info) + call field2_desc%clone(grid_desc%descs(1,2), info) + call field1_desc%clone(grid_desc%descs(2,1), info) + call field2_desc%clone(grid_desc%descs(2,2), info) + + !--------------------------------------------------------------- + ! 7) composed global operator (what CG will use as its matrix) + !--------------------------------------------------------------- + call psb_cd_nest_compose(grid_desc, desc_global, info) + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: psb_cd_nest_compose info=', info + goto 9999 + end if + call psb_d_nest_base_setup(nest_operator, block_storage, grid_desc, desc_global, info) + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: psb_d_nest_base_setup info=', info + goto 9999 + end if + allocate(global_operator%a, source=nest_operator) + call global_operator%set_nrows(desc_global%get_local_rows()) + call global_operator%set_ncols(desc_global%get_local_cols()) + call global_operator%set_asb() + + !--------------------------------------------------------------- + ! 8) consistent RHS: x_exact = 1, rhs = M * x_exact (via the nested operator) + !--------------------------------------------------------------- + call psb_geall(x_exact, desc_global, info) + do i_local_row = 1, desc_global%get_local_rows() + call desc_global%l2g(i_local_row, field1_global_row, info) + insert_value(1) = 1.0_psb_dpk_ + call psb_geins(1, [field1_global_row], insert_value, x_exact, desc_global, info) + end do + call psb_geasb(x_exact, desc_global, info) + + call psb_geall(rhs, desc_global, info); call psb_geasb(rhs, desc_global, info) + call psb_spmm(done, global_operator, x_exact, dzero, rhs, desc_global, info) + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (RHS) info=', info + goto 9999 + end if + + norm_x_exact = psb_genrm2(x_exact, desc_global, info) + + !--------------------------------------------------------------- + ! 9) identity preconditioner (NONE): CG exercises only the operator + !--------------------------------------------------------------- + call preconditioner%init(context, 'NONE', info) + call preconditioner%build(global_operator, desc_global, info) + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: preconditioner%build info=', info + goto 9999 + end if + + !--------------------------------------------------------------- + ! 10) solve with the standard PSBLAS CG + !--------------------------------------------------------------- + call psb_geall(x_solution, desc_global, info); call psb_geasb(x_solution, desc_global, info) + call psb_krylov('CG', global_operator, preconditioner, rhs, x_solution, stop_tol, desc_global, info, & + & itmax=max_iter, iter=n_iter, err=final_residual, itrace=trace_level, istop=stop_criterion) + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG) info=', info + goto 9999 + end if + + !--------------------------------------------------------------- + ! 11) solution error: || x_solution - x_exact || / || x_exact || + !--------------------------------------------------------------- + call psb_geaxpby(-done, x_exact, done, x_solution, desc_global, info) ! x_solution <- x_solution - x_exact + solution_error = psb_genrm2(x_solution, desc_global, info) / norm_x_exact + + if (my_rank == 0) then + write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(global)=', 2*field_size + write(*,'(a,i0)') ' CG iterations = ', n_iter + write(*,'(a,es12.4)') ' CG relative residual = ', final_residual + write(*,'(a,es12.4)') ' ||x - x_exact||/||x_ex|| = ', solution_error + if ((n_iter < max_iter) .and. (solution_error <= solution_tol)) then + write(*,*) '[PASS] CG converges on the global nested operator' + else + write(*,*) '[FAIL] CG does not converge / wrong solution (tol ', solution_tol, ')' + end if + end if + +9999 continue + call psb_exit(context) + +end program psb_d_nest_cg_test diff --git a/test/nested/psb_d_nest_glob_test.F90 b/test/nested/psb_d_nest_glob_test.F90 new file mode 100644 index 000000000..80054d09e --- /dev/null +++ b/test/nested/psb_d_nest_glob_test.F90 @@ -0,0 +1,218 @@ +! +! 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 prior 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_nest_glob_test.F90 +! +! Program: psb_d_nest_glob_test +! Author: Simone Staccone (Stack-1) +! +! Validates the "global nested operator" path built through the +! psb_d_nest_matrix utility (init/ins/asb): the user only supplies the field +! sizes and the block values, and obtains nested_matrix%a_glob / +! nested_matrix%desc_glob ready for psb_spmm. The result is compared against +! the SAME matrix assembled monolithically in CSR on the same global +! descriptor (oracle). +! +! 2x2 operator (fields of size field_size): +! [ A B^T ] A = tridiag(-1, 2, -1) (block 1,1) +! [ B 0 ] B^T = 0.5 * I (block 1,2) +! B = 0.3 * I (block 2,1) +! (2,2) absent +! +! Run: ./psb_d_nest_glob_test (serial) +! mpirun -np 4 ./psb_d_nest_glob_test +! +program psb_d_nest_glob_test + use psb_base_mod + use psb_util_mod + use psb_d_nest_mod + implicit none + + type(psb_ctxt_type) :: context + integer(psb_ipk_) :: my_rank, num_procs, info, i_local_row + integer(psb_ipk_) :: entry_idx, field1_local_rows, field2_local_rows + integer(psb_lpk_) :: global_row, global_col, field_size + + type(psb_d_nest_matrix) :: nested_matrix ! the nested operator (init/ins/asb) + type(psb_dspmat_type) :: monolithic_ref ! monolithic CSR oracle + type(psb_d_vect_type) :: x_vec, y_nested, y_monolithic + + integer(psb_lpk_), allocatable :: entry_rows(:), entry_cols(:) + real(psb_dpk_), allocatable :: entry_vals(:) + real(psb_dpk_) :: insert_value(1) + real(psb_dpk_) :: mismatch_norm + real(psb_dpk_), parameter :: tolerance = 1.0e-10_psb_dpk_ + + call psb_init(context) + call psb_info(context, my_rank, num_procs) + + field_size = 32 ! global size of each field + + !--------------------------------------------------------------- + ! 1) build the 2x2 nested operator through the utility + !--------------------------------------------------------------- + call nested_matrix%init(context, [field_size, field_size], info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL: nested_matrix%init info=', info; goto 9999 + end if + field1_local_rows = nested_matrix%field_desc(1)%get_local_rows() + field2_local_rows = nested_matrix%field_desc(2)%get_local_rows() + + !--------------------------------------------------------------- + ! 2) insert the block values (owned rows only) + !--------------------------------------------------------------- + ! A = tridiag(-1,2,-1) -> block (1,1) + allocate(entry_rows(3*field1_local_rows), entry_cols(3*field1_local_rows), & + & entry_vals(3*field1_local_rows)) + entry_idx = 0 + do i_local_row = 1, field1_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = global_row + entry_cols(entry_idx) = global_row + entry_vals(entry_idx) = 2.0_psb_dpk_ + if (global_row > 1) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = global_row + entry_cols(entry_idx) = global_row - 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + if (global_row < field_size) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = global_row + entry_cols(entry_idx) = global_row + 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call nested_matrix%ins(1, 1, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! B^T = 0.5 I -> block (1,2): rows in field 1, columns in field 2 + allocate(entry_rows(field1_local_rows), entry_cols(field1_local_rows), entry_vals(field1_local_rows)) + entry_idx = 0 + do i_local_row = 1, field1_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = global_row + entry_cols(entry_idx) = global_row + entry_vals(entry_idx) = 0.5_psb_dpk_ + end do + call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! B = 0.3 I -> block (2,1): rows in field 2, columns in field 1 + allocate(entry_rows(field2_local_rows), entry_cols(field2_local_rows), entry_vals(field2_local_rows)) + entry_idx = 0 + do i_local_row = 1, field2_local_rows + call nested_matrix%field_desc(2)%l2g(i_local_row, global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = global_row + entry_cols(entry_idx) = global_row + entry_vals(entry_idx) = 0.3_psb_dpk_ + end do + call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + call nested_matrix%asb(info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL: nested_matrix%asb info=', info; goto 9999 + end if + + !--------------------------------------------------------------- + ! 3) monolithic oracle on nested_matrix%desc_glob (global offsets: + ! field 1 -> g ; field 2 -> field_size + g) + !--------------------------------------------------------------- + call psb_spall(monolithic_ref, nested_matrix%desc_glob, info, & + & nnz=5*nested_matrix%desc_glob%get_local_rows()) + do i_local_row = 1, field1_local_rows ! field-1 rows + call nested_matrix%field_desc(1)%l2g(i_local_row, global_row, info) + insert_value(1) = 2.0_psb_dpk_ + call psb_spins(1,[global_row],[global_row],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + if (global_row > 1) then + insert_value(1)=-1.0_psb_dpk_ + call psb_spins(1,[global_row],[global_row-1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + end if + if (global_row < field_size) then + insert_value(1)=-1.0_psb_dpk_ + call psb_spins(1,[global_row],[global_row+1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + end if + global_col = field_size + global_row + insert_value(1) = 0.5_psb_dpk_ ! B^T + call psb_spins(1,[global_row],[global_col],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + end do + do i_local_row = 1, field2_local_rows ! field-2 rows + call nested_matrix%field_desc(2)%l2g(i_local_row, global_row, info) + global_col = global_row + insert_value(1) = 0.3_psb_dpk_ ! B + call psb_spins(1,[field_size+global_row],[global_col],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + end do + call psb_spasb(monolithic_ref, nested_matrix%desc_glob, info, dupl=psb_dupl_add_) + + !--------------------------------------------------------------- + ! 4) compare the two matrix-vector products on a distinct-valued x (x[g] = g) + !--------------------------------------------------------------- + call psb_geall(x_vec, nested_matrix%desc_glob, info) + do i_local_row = 1, nested_matrix%desc_glob%get_local_rows() + call nested_matrix%desc_glob%l2g(i_local_row, global_row, info) + insert_value(1) = real(global_row, psb_dpk_) + call psb_geins(1, [global_row], insert_value, x_vec, nested_matrix%desc_glob, info) + end do + call psb_geasb(x_vec, nested_matrix%desc_glob, info) + + call psb_geall(y_nested, nested_matrix%desc_glob, info); call psb_geasb(y_nested, nested_matrix%desc_glob, info) + call psb_geall(y_monolithic, nested_matrix%desc_glob, info); call psb_geasb(y_monolithic, nested_matrix%desc_glob, info) + + call psb_spmm(done, nested_matrix%a_glob, x_vec, dzero, y_nested, nested_matrix%desc_glob, info) ! via nested csmv + if (info /= psb_success_) then + if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (nested) info=', info + goto 9999 + end if + call psb_spmm(done, monolithic_ref, x_vec, dzero, y_monolithic, nested_matrix%desc_glob, info) ! CSR oracle + + call psb_geaxpby(done, y_nested, -done, y_monolithic, nested_matrix%desc_glob, info) + mismatch_norm = psb_genrm2(y_monolithic, nested_matrix%desc_glob, info) + + if (my_rank == 0) then + write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(field)=', field_size + write(*,'(a,es12.4)') ' ||nested - monolithic||_2 = ', mismatch_norm + if (mismatch_norm <= tolerance) then + write(*,*) '[PASS] nested global operator matches monolithic CSR' + else + write(*,*) '[FAIL] mismatch above tolerance ', tolerance + end if + end if + + call nested_matrix%free(info) + +9999 continue + call psb_exit(context) + +end program psb_d_nest_glob_test diff --git a/test/nested/psb_d_nest_rect_test.F90 b/test/nested/psb_d_nest_rect_test.F90 new file mode 100644 index 000000000..af88d2566 --- /dev/null +++ b/test/nested/psb_d_nest_rect_test.F90 @@ -0,0 +1,217 @@ +! +! 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 prior 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_nest_rect_test.F90 +! +! Program: psb_d_nest_rect_test +! Author: Simone Staccone (Stack-1) +! +! Like psb_d_nest_glob_test but with fields of DIFFERENT size (|V| = 2|Q|) and +! GENUINELY RECTANGULAR off-diagonal blocks. The operator is built with the +! psb_d_nest_matrix utility (init/ins/asb): the user only inserts the values, +! while the cross-field halo registration and the construction of the +! rectangular local blocks are handled internally. Compared against the +! monolithic CSR oracle. +! +! [ A B^T ] A : V x V tridiag(-1,2,-1) +! [ B 0 ] B^T : V x Q rectangular (row r -> col mod(r-1,nQ)+1, val 0.5) +! B : Q x V rectangular (row q -> cols q and q+nQ, val 0.3) +! (2,2) absent +! +! Run: ./psb_d_nest_rect_test ; mpirun -np 4 ./psb_d_nest_rect_test +! +program psb_d_nest_rect_test + use psb_base_mod + use psb_util_mod + use psb_d_nest_mod + implicit none + + type(psb_ctxt_type) :: context + integer(psb_ipk_) :: my_rank, num_procs, info, i_local_row + integer(psb_ipk_) :: entry_idx, v_local_rows, q_local_rows + integer(psb_lpk_) :: v_global_row, q_global_row, q_col, v_size, q_size + + type(psb_d_nest_matrix) :: nested_matrix + type(psb_dspmat_type) :: monolithic_ref + type(psb_d_vect_type) :: x_vec, y_nested, y_monolithic + real(psb_dpk_) :: insert_value(1) + + integer(psb_lpk_), allocatable :: entry_rows(:), entry_cols(:) + real(psb_dpk_), allocatable :: entry_vals(:) + real(psb_dpk_) :: mismatch_norm + real(psb_dpk_), parameter :: tolerance = 1.0e-10_psb_dpk_ + + call psb_init(context) + call psb_info(context, my_rank, num_procs) + + q_size = 16 ! global size of field Q + v_size = 2*q_size ! global size of field V (|V| = 2|Q|) + + !--------------------------------------------------------------- + ! 1) build the 2x2 nested operator (fields V, Q) + !--------------------------------------------------------------- + call nested_matrix%init(context, [v_size, q_size], info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL: nested_matrix%init info=', info; goto 9999 + end if + v_local_rows = nested_matrix%field_desc(1)%get_local_rows() + q_local_rows = nested_matrix%field_desc(2)%get_local_rows() + + !--------------------------------------------------------------- + ! 2) insert the blocks (owned rows only) + !--------------------------------------------------------------- + ! A = tridiag(-1,2,-1) -> block (1,1), V x V + allocate(entry_rows(3*v_local_rows), entry_cols(3*v_local_rows), entry_vals(3*v_local_rows)) + entry_idx = 0 + do i_local_row = 1, v_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, v_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = v_global_row + entry_cols(entry_idx) = v_global_row + entry_vals(entry_idx) = 2.0_psb_dpk_ + if (v_global_row > 1) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = v_global_row + entry_cols(entry_idx) = v_global_row - 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + if (v_global_row < v_size) then + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = v_global_row + entry_cols(entry_idx) = v_global_row + 1_psb_lpk_ + entry_vals(entry_idx) = -1.0_psb_dpk_ + end if + end do + call nested_matrix%ins(1, 1, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! B^T rectangular -> block (1,2), V x Q : row r -> col mod(r-1,nQ)+1, val 0.5 + allocate(entry_rows(v_local_rows), entry_cols(v_local_rows), entry_vals(v_local_rows)) + entry_idx = 0 + do i_local_row = 1, v_local_rows + call nested_matrix%field_desc(1)%l2g(i_local_row, v_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = v_global_row + entry_cols(entry_idx) = mod(v_global_row-1_psb_lpk_, q_size)+1 + entry_vals(entry_idx) = 0.5_psb_dpk_ + end do + call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + ! B rectangular -> block (2,1), Q x V : row q -> cols q and q+nQ, val 0.3 + allocate(entry_rows(2*q_local_rows), entry_cols(2*q_local_rows), entry_vals(2*q_local_rows)) + entry_idx = 0 + do i_local_row = 1, q_local_rows + call nested_matrix%field_desc(2)%l2g(i_local_row, q_global_row, info) + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = q_global_row + entry_cols(entry_idx) = q_global_row + entry_vals(entry_idx) = 0.3_psb_dpk_ + entry_idx = entry_idx + 1 + entry_rows(entry_idx) = q_global_row + entry_cols(entry_idx) = q_global_row + q_size + entry_vals(entry_idx) = 0.3_psb_dpk_ + end do + call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info) + deallocate(entry_rows, entry_cols, entry_vals) + + call nested_matrix%asb(info) + if (info /= psb_success_) then + if (my_rank==0) write(*,*) 'FAIL: nested_matrix%asb info=', info; goto 9999 + end if + + !--------------------------------------------------------------- + ! 3) monolithic oracle on nested_matrix%desc_glob (global offsets: V -> g, Q -> v_size + g) + !--------------------------------------------------------------- + call psb_spall(monolithic_ref, nested_matrix%desc_glob, info, & + & nnz=6*nested_matrix%desc_glob%get_local_rows()) + do i_local_row = 1, v_local_rows ! V rows + call nested_matrix%field_desc(1)%l2g(i_local_row, v_global_row, info) + insert_value(1)=2.0_psb_dpk_ + call psb_spins(1,[v_global_row],[v_global_row],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + if (v_global_row>1) then + insert_value(1)=-1.0_psb_dpk_ + call psb_spins(1,[v_global_row],[v_global_row-1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info) + end if + if (v_global_row=c, 0 otherwise -! T16 psb_d_nest_mask mask operation; t=.true. if all satisfied -! T17 psb_d_nest_upd_xyz y=alpha*x+beta*y; z=gamma*y_new+delta*z -! -program psb_d_pde_nest_psblas - use psb_base_mod - use psb_desc_nest_mod - use psb_d_nest_mod - use psb_d_nest_tools_mod, only : psb_geall_nest, psb_geasb_nest, psb_gefree_nest, & - & psb_geins_nest, psb_spall_nest, psb_spins_nest, & - & psb_spasb_nest - implicit none - - !------------------------------------------------------------------ - ! Parallel context - !------------------------------------------------------------------ - type(psb_ctxt_type) :: ctxt - integer(psb_ipk_) :: iam, np, info - - !------------------------------------------------------------------ - ! Problem size - !------------------------------------------------------------------ - integer(psb_ipk_), parameter :: n = 100 - integer(psb_ipk_) :: nlr, iloc, i - integer(psb_lpk_) :: iglob - - !------------------------------------------------------------------ - ! Per-block descriptors (identical layout as psb_d_pde_nest_first) - !------------------------------------------------------------------ - type(psb_desc_type) :: desc1, desc2, desc3, desc4 - - !------------------------------------------------------------------ - ! Nested descriptor and nested sparse matrix - !------------------------------------------------------------------ - type(psb_desc_nest_type) :: descs - type(psb_d_nest_sparse_mat) :: anest - - !------------------------------------------------------------------ - ! Individual sparse matrices (A11 = Laplacian, A12 = I, A21 = I) - !------------------------------------------------------------------ - type(psb_dspmat_type) :: a11, a12, a21 - - !------------------------------------------------------------------ - ! Work nested vectors (xnest, ynest, znest reused across tests) - !------------------------------------------------------------------ - type(psb_d_nest_vect_type) :: xnest, ynest, znest - - !------------------------------------------------------------------ - ! Insert buffers - !------------------------------------------------------------------ - integer(psb_lpk_) :: grow(1) - real(psb_dpk_) :: gval(1) - - !------------------------------------------------------------------ - ! Scalar results and expected values - !------------------------------------------------------------------ - real(psb_dpk_) :: res, res2, expected - real(psb_dpk_), parameter :: tol = 1.0e-10_psb_dpk_ - logical :: t_mask - - !------------------------------------------------------------------ - ! Test pass/fail counter - !------------------------------------------------------------------ - integer(psb_ipk_) :: npass, nfail - - character(len=40) :: name = 'psb_d_pde_nest_psblas' - - !================================================================== - ! Initialise MPI - !================================================================== - call psb_init(ctxt) - call psb_info(ctxt, iam, np) - call psb_set_errverbosity(itwo) - - npass = 0 - nfail = 0 - - !================================================================== - ! Build per-block descriptors - ! Use exact block distribution so total rows = n exactly. - ! Ceiling division (n+np-1)/np gives nlr*np > n phantom rows - ! when np does not divide n evenly. - !================================================================== - nlr = n / np - if (iam < mod(n, np)) nlr = nlr + 1_psb_ipk_ - nlr = max(1_psb_ipk_, nlr) - call psb_cdall(ctxt, desc1, info, nl=nlr) - call psb_cdall(ctxt, desc2, info, nl=nlr) - call desc1%clone(desc3, info) - - !================================================================== - ! A(1,1) = tridiagonal Laplacian - !================================================================== - call psb_spall(a11, desc1, info, nnz=3*desc1%get_local_rows()) - do iloc = 1, desc1%get_local_rows() - call desc1%l2g(iloc, iglob, info) - grow(1)=iglob; gval(1)=2.0_psb_dpk_ - call psb_spins(1,grow,grow,gval,a11,desc1,info) - if (iglob>1) then - grow(1)=iglob; gval(1)=-1.0_psb_dpk_ - call psb_spins(1,grow,[iglob-1_psb_lpk_],gval,a11,desc1,info) - end if - if (iglob y1=2 (2 boundary rows) - ! rows 2..n-1 give Lap=0, I=1 => y1=1 (n-2 interior rows) - ! Block 2: I*ones = 1 for all n rows - ! amax=2, gemin=1, geasum = (n+2) + n = 2n+2 - res = psb_d_nest_geamax(ynest, descs, info) - call check('T01 spmm amax(y)=2', res, 2.0_psb_dpk_, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T01 spmm gemin(y)=1', res, done, tol, npass, nfail, iam) - res = psb_d_nest_geasum(ynest, descs, info) - expected = 2.0_psb_dpk_ * real(n, psb_dpk_) + 2.0_psb_dpk_ - call check('T01 spmm geasum(y)=2n+2', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T02: psb_d_nest_geaxpby - ! x = all 3s, y = all 2s - ! y = 2*x + (-1)*y => y = 2*3 - 2 = 4 (all 4s) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T02: psb_d_nest_geaxpby (y = 2*x - y, x=3 y=2 => 4)' - - call set_nest_val(xnest, 3.0_psb_dpk_) - call set_nest_val(ynest, 2.0_psb_dpk_) - - call psb_d_nest_geaxpby(2.0_psb_dpk_, xnest, -done, ynest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 4.0_psb_dpk_ - res = psb_d_nest_geamax(ynest, descs, info) - call check('T02 geaxpby amax(y)=4', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T02 geaxpby amin(y)=4', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T03: psb_d_nest_genrm2 - ! x = all 1s => ||x||_2 = sqrt(2*n) = sqrt(20) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T03: psb_d_nest_genrm2 (x=1 => sqrt(2n))' - - call set_nest_val(xnest, done) - - res = psb_d_nest_genrm2(xnest, descs, info) - expected = sqrt(2.0_psb_dpk_ * real(n, psb_dpk_)) - call check('T03 genrm2(ones)', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T04: psb_d_nest_genrm2s (subroutine form; result must equal T03) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T04: psb_d_nest_genrm2s (subroutine form of genrm2)' - - call psb_d_nest_genrm2s(res2, xnest, descs, info) - call check('T04 genrm2s == genrm2', res2, res, tol, npass, nfail, iam) - - !================================================================== - ! T05: psb_d_nest_gedot - ! x = all 1s, y = all 2s => dot = 2 * 2*n = 40 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T05: psb_d_nest_gedot (x=1 y=2 => 2*2n=40)' - - call set_nest_val(xnest, done) - call set_nest_val(ynest, 2.0_psb_dpk_) - - res = psb_d_nest_gedot(xnest, ynest, descs, info) - expected = 2.0_psb_dpk_ * 2.0_psb_dpk_ * real(n, psb_dpk_) - call check('T05 gedot', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T06: psb_d_nest_geamax - ! x = all 5s => ||x||_inf = 5 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T06: psb_d_nest_geamax (x=5 => 5)' - - call set_nest_val(xnest, 5.0_psb_dpk_) - - res = psb_d_nest_geamax(xnest, descs, info) - expected = 5.0_psb_dpk_ - call check('T06 geamax', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T07: psb_d_nest_geasum - ! x = all 1s => ||x||_1 = 2*n = 20 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T07: psb_d_nest_geasum (x=1 => 2n=20)' - - call set_nest_val(xnest, done) - - res = psb_d_nest_geasum(xnest, descs, info) - expected = 2.0_psb_dpk_ * real(n, psb_dpk_) - call check('T07 geasum', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T08: psb_d_nest_gemin - ! x = all 7s => min = 7 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T08: psb_d_nest_gemin (x=7 => 7)' - - call set_nest_val(xnest, 7.0_psb_dpk_) - - res = psb_d_nest_gemin(xnest, descs, info) - expected = 7.0_psb_dpk_ - call check('T08 gemin', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T09: psb_d_nest_minquotient - ! x = all 3s, y = all 6s => min(x/y) = 0.5 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T09: psb_d_nest_minquotient (x=3 y=6 => 0.5)' - - call set_nest_val(xnest, 3.0_psb_dpk_) - call set_nest_val(ynest, 6.0_psb_dpk_) - - res = psb_d_nest_minquotient(xnest, ynest, descs, info) - expected = 0.5_psb_dpk_ - call check('T09 minquotient', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T10: psb_d_nest_gemlt - ! x = all 2s, y = all 4s => y = y * x = 8 (result in y) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T10: psb_d_nest_gemlt (x=2 y=4 => y=8)' - - call set_nest_val(xnest, 2.0_psb_dpk_) - call set_nest_val(ynest, 4.0_psb_dpk_) - - call psb_d_nest_gemlt(xnest, ynest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 8.0_psb_dpk_ - res = psb_d_nest_geamax(ynest, descs, info) - call check('T10 gemlt amax(y)=8', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T10 gemlt amin(y)=8', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T11: psb_d_nest_gediv - ! x = all 6s, y = all 3s => x = x / y = 2 (result in x) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T11: psb_d_nest_gediv (x=6 y=3 => x=2)' - - call set_nest_val(xnest, 6.0_psb_dpk_) - call set_nest_val(ynest, 3.0_psb_dpk_) - - call psb_d_nest_gediv(xnest, ynest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 2.0_psb_dpk_ - res = psb_d_nest_geamax(xnest, descs, info) - call check('T11 gediv amax(x)=2', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(xnest, descs, info) - call check('T11 gediv amin(x)=2', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T12: psb_d_nest_geinv - ! x = all 4s => y = 1/x = 0.25 (result in y) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T12: psb_d_nest_geinv (x=4 => y=0.25)' - - call set_nest_val(xnest, 4.0_psb_dpk_) - - call psb_d_nest_geinv(xnest, ynest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 0.25_psb_dpk_ - res = psb_d_nest_geamax(ynest, descs, info) - call check('T12 geinv amax(y)=0.25', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T12 geinv amin(y)=0.25', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T13: psb_d_nest_geabs - ! x = all -3s => y = |x| = 3 (result in y) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T13: psb_d_nest_geabs (x=-3 => y=3)' - - call set_nest_val(xnest, -3.0_psb_dpk_) - - call psb_d_nest_geabs(xnest, ynest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 3.0_psb_dpk_ - res = psb_d_nest_geamax(ynest, descs, info) - call check('T13 geabs amax(y)=3', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T13 geabs amin(y)=3', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T14: psb_d_nest_geaddconst - ! x = all 2s, b = 7.0 => z = x + 7 = 9 (result in z) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T14: psb_d_nest_geaddconst (x=2 b=7 => z=9)' - - call set_nest_val(xnest, 2.0_psb_dpk_) - - call psb_d_nest_geaddconst(xnest, 7.0_psb_dpk_, znest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 9.0_psb_dpk_ - res = psb_d_nest_geamax(znest, descs, info) - call check('T14 geaddconst amax(z)=9', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(znest, descs, info) - call check('T14 geaddconst amin(z)=9', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T15a: psb_d_nest_gecmp — entries satisfy threshold - ! x = all 3s, c = 2.0 => z(i)=1 (since |3| >= 2) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T15a: psb_d_nest_gecmp (x=3 c=2 => z=1)' - - call set_nest_val(xnest, 3.0_psb_dpk_) - - call psb_d_nest_gecmp(xnest, 2.0_psb_dpk_, znest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = done - res = psb_d_nest_geamax(znest, descs, info) - call check('T15a gecmp amax(z)=1', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(znest, descs, info) - call check('T15a gecmp amin(z)=1', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T15b: psb_d_nest_gecmp — entries do not satisfy threshold - ! x = all 1s, c = 2.0 => z(i)=0 (since |1| < 2) - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T15b: psb_d_nest_gecmp (x=1 c=2 => z=0)' - - call set_nest_val(xnest, done) - - call psb_d_nest_gecmp(xnest, 2.0_psb_dpk_, znest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = dzero - res = psb_d_nest_geamax(znest, descs, info) - call check('T15b gecmp amax(z)=0', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! T16: psb_d_nest_mask - ! Semantics: mask(c, x, m, t) - ! c = values to test (first arg) - ! x = constraint-type indicators (second arg): - ! 2 => satisfied if c(i) > 0 - ! 1 => satisfied if c(i) >= 0 - ! -1 => satisfied if c(i) <= 0 - ! -2 => satisfied if c(i) < 0 - ! m = output mask (0=satisfied, 1=violated) - ! t = .true. iff all entries satisfied - ! - ! Case: c = all +3 (positive), x = all 2 (check > 0) => m=0 t=T - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') 'T16: psb_d_nest_mask (c=3 x=2 => m=0 t=.true.)' - - call set_nest_val(xnest, 3.0_psb_dpk_) ! values to test - call set_nest_val(ynest, 2.0_psb_dpk_) ! constraint indicators (type 2: check > 0) - - call psb_d_nest_mask(xnest, ynest, znest, t_mask, descs, info) - if (info /= psb_success_) goto 9999 - - if (iam == 0) then - if (t_mask) then - write(*,'(a)') ' T16 mask: t=.true. PASS (all constraints satisfied)' - npass = npass + 1 - else - write(*,'(a)') ' T16 mask: t=.false. FAIL (expected .true.)' - nfail = nfail + 1 - end if - end if - - !------------------------------------------------------------------ - ! T16b: c = all -3 (negative), x = all 2 (check > 0) => m=1 t=F - !------------------------------------------------------------------ - if (iam == 0) write(*,'(a)') 'T16b: psb_d_nest_mask (c=-3 x=2 => m=1 t=.false.)' - - call set_nest_val(xnest, -3.0_psb_dpk_) ! values (negative) - call set_nest_val(ynest, 2.0_psb_dpk_) ! indicators (type 2: check > 0) - - call psb_d_nest_mask(xnest, ynest, znest, t_mask, descs, info) - if (info /= psb_success_) goto 9999 - - if (iam == 0) then - if (.not. t_mask) then - write(*,'(a)') ' T16b mask: t=.false. PASS (all constraints violated)' - npass = npass + 1 - else - write(*,'(a)') ' T16b mask: t=.true. FAIL (expected .false.)' - nfail = nfail + 1 - end if - end if - - !================================================================== - ! T17: psb_d_nest_upd_xyz - ! Computes: y_new = alpha*x + beta*y - ! z_new = gamma*y_new + delta*z - ! - ! x=1, y=2, z=3, alpha=2, beta=3, gamma=4, delta=5 - ! => y_new = 2*1 + 3*2 = 8 - ! => z_new = 4*8 + 5*3 = 47 - !================================================================== - if (iam == 0) write(*,'(/,a)') repeat('=',60) - if (iam == 0) write(*,'(a)') & - 'T17: psb_d_nest_upd_xyz (x=1 y=2 z=3 a=2 b=3 g=4 d=5 => y=8 z=47)' - - call set_nest_val(xnest, done) - call set_nest_val(ynest, 2.0_psb_dpk_) - call set_nest_val(znest, 3.0_psb_dpk_) - - call psb_d_nest_upd_xyz(2.0_psb_dpk_, 3.0_psb_dpk_, & - & 4.0_psb_dpk_, 5.0_psb_dpk_, & - & xnest, ynest, znest, descs, info) - if (info /= psb_success_) goto 9999 - - expected = 8.0_psb_dpk_ - res = psb_d_nest_geamax(ynest, descs, info) - call check('T17 upd_xyz amax(y)=8', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(ynest, descs, info) - call check('T17 upd_xyz amin(y)=8', res, expected, tol, npass, nfail, iam) - - expected = 47.0_psb_dpk_ - res = psb_d_nest_geamax(znest, descs, info) - call check('T17 upd_xyz amax(z)=47', res, expected, tol, npass, nfail, iam) - res = psb_d_nest_gemin(znest, descs, info) - call check('T17 upd_xyz amin(z)=47', res, expected, tol, npass, nfail, iam) - - !================================================================== - ! Summary - !================================================================== - if (iam == 0) then - write(*,'(/,a)') repeat('=',60) - write(*,'(a,i0,a,i0,a)') & - ' RESULTS: ', npass, ' passed, ', nfail, ' failed' - write(*,'(a)') repeat('=',60) - end if - - !================================================================== - ! Clean up - !================================================================== - call psb_gefree_nest(xnest, descs, info) - call psb_gefree_nest(ynest, descs, info) - call psb_gefree_nest(znest, descs, info) - - call psb_cdfree(desc1, info) - call psb_cdfree(desc2, info) - call psb_cdfree(desc3, info) - call psb_cdfree(desc4, info) - - call psb_exit(ctxt) - stop - -9999 continue - write(psb_err_unit,*) trim(name), ': error info=', info, ' rank=', iam - call psb_error(ctxt) - call psb_exit(ctxt) - stop - -contains - - !------------------------------------------------------------------ - ! Set every local entry of every block to val - !------------------------------------------------------------------ - subroutine set_nest_val(v, val) - use psb_base_mod - type(psb_d_nest_vect_type), intent(inout) :: v - real(psb_dpk_), intent(in) :: val - integer(psb_ipk_) :: k, linfo - linfo = 0 - do k = 1, v%nblocks - call v%vects(k)%set(val, linfo) - end do - end subroutine set_nest_val - - !------------------------------------------------------------------ - ! Scalar pass/fail check with tolerance - !------------------------------------------------------------------ - subroutine check(label, got, expected, tol, np_, nf_, myrank) - use psb_base_mod - character(len=*), intent(in) :: label - real(psb_dpk_), intent(in) :: got, expected, tol - integer(psb_ipk_), intent(inout) :: np_, nf_ - integer(psb_ipk_), intent(in) :: myrank - - if (myrank /= 0) return - if (abs(got - expected) <= tol * max(done, abs(expected))) then - write(*,'(2x,a,a,f16.10,a,f16.10)') & - 'PASS ', trim(label)//' got=', got, ' exp=', expected - np_ = np_ + 1 - else - write(*,'(2x,a,a,f16.10,a,f16.10)') & - 'FAIL ', trim(label)//' got=', got, ' exp=', expected - nf_ = nf_ + 1 - end if - end subroutine check - - !------------------------------------------------------------------ - ! Print every block of a nested vector (one rank at a time). - ! Each process flushes stdout before the barrier so that buffered - ! output does not bleed into the next process's print window. - !------------------------------------------------------------------ - subroutine print_nest_vec(v, label, myrank, nprocs, myctxt, ds) - use psb_base_mod - use iso_fortran_env, only: output_unit - type(psb_d_nest_vect_type), intent(inout) :: v - character(len=*), intent(in) :: label - integer(psb_ipk_), intent(in) :: myrank, nprocs - type(psb_ctxt_type), intent(in) :: myctxt - type(psb_desc_nest_type), intent(in) :: ds - - integer(psb_ipk_) :: blk, ip, k, nr, linfo - real(psb_dpk_), allocatable :: vals(:) - - do blk = 1, v%nblocks - nr = ds%descs(blk,blk)%get_local_rows() - do ip = 0, nprocs-1 - call psb_barrier(myctxt) - if (myrank == ip) then - write(*,'(a,a,a,i0,a)') ' [', trim(label), '] block ', blk, ':' - linfo = 0 - allocate(vals(nr), stat=linfo) - if (linfo == 0) vals = v%vects(blk)%get_vect() - do k = 1, nr - write(*,'(4x,i4,f14.6)') k, vals(k) - end do - deallocate(vals) - flush(output_unit) - end if - end do - call psb_barrier(myctxt) - end do - end subroutine print_nest_vec - -end program psb_d_pde_nest_psblas