fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-20 21:02:26 +00:00
Code Issues Packages Projects Releases Wiki Activity

[ADD] Nested (MATNEST-style) block matrix support and tests

Browse Source 
Add a block-structured distributed operator that presents itself to Krylov
solvers and preconditioners as a single ordinary distributed matrix (the
PSBLAS analogue of PETSc MATNEST), targeting saddle-point systems
M = [[A, B^T], [B, 0]] with possibly rectangular sub-blocks.

Library (base/modules):
- psb_desc_nest_mod, psb_d_nest_mat_mod: grid of per-field descriptors and
  per-block sparse storage.
- psb_d_nest_base_mat_mod: psb_d_nest_base_mat, the operator extending
  psb_d_base_sparse_mat (local csmv, free, field-split hooks for a future
  block preconditioner).
- psb_cd_nest_tools_mod / psb_d_nest_tools_mod: composed global descriptor
  with union halo (psb_cd_nest_compose) and rectangular local block builder
  (psb_d_nest_rect_block), plus the per-block assembly wrappers.
- psb_d_nest_builder_mod: psb_d_nest_matrix, the user frontend with the
  init/ins/asb/free pattern hiding all descriptor/halo/compose/setup
  boilerplate.
- psb_d_nest_mod: umbrella module (use psb_d_nest_mod).

Remove the earlier bespoke per-block prototype (comm/psblas/vect modules and
the pde_nest_psblas test) superseded by the single MATNEST design.

Tests (test/nested): glob (square operator vs monolithic CSR oracle), rect
(genuinely rectangular blocks), cg (low-level path, ill-conditioned SPD
red-black Laplacian solved with standard CG), builder (same solve via the
utility), plus a README describing the design and usage. All pass serially
and in parallel, with results invariant to the process count.

Build hooks updated (autotools Makefiles + CMakeLists); the nested tests are
relocated out of test/pdegen into test/nested.

Author: Simone Staccone (Stack-1)
This commit is contained in:
Stack-1
2026-06-11 16:53:46 +02:00
parent 784c3cc0b4
commit acdd2e9eb5
23 changed files with 2384 additions and 1889 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/Makefile
+14 -19
View File
@@ -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)
base/modules/comm/psb_d_nest_comm_mod.f90
-93
View File
@@ -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
base/modules/desc/psb_desc_nest_mod.f90
+43 -23
View File
@@ -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
base/modules/psb_d_nest_mod.f90
+7 -4
View File
@@ -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
base/modules/psblas/psb_d_nest_psblas_mod.f90
-626
View File
@@ -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
base/modules/serial/psb_d_nest_base_mat_mod.F90
+405
View File
@@ -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
base/modules/serial/psb_d_nest_mat_mod.f90
+31 -31
View File
@@ -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
base/modules/serial/psb_d_nest_vect_mod.f90
-109
View File
@@ -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
base/modules/tools/psb_cd_nest_tools_mod.F90
+202 -75
View File
@@ -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<k} n_m
!
! Each process owns its slice of every field; the global halo is the union of
! the per-field halos, each remapped by its field offset. Once composed, the
! nested operator can be presented to Krylov/AMG4PSBLAS as a standard
! distributed matrix/vector (MATNEST-style).
!
! Assumes a square block structure (nrblocks == ncblocks); field k is taken to
! be column k, whose distribution and halo are read from descs(1,k) (all
! descs(i,k) for fixed k share the same column space).
!
subroutine psb_cd_nest_compose(desc_grid, desc_global, info)
type(psb_desc_nest_type), intent(in) :: desc_grid
type(psb_desc_type), intent(out) :: desc_global
integer(psb_ipk_), intent(out) :: info
type(psb_ctxt_type) :: ctxt
integer(psb_ipk_) :: n_fields, i_field, i_loc, n_owned, n_local, owned_count, halo_count
integer(psb_lpk_) :: global_idx
integer(psb_lpk_), allocatable :: field_offset(:), owned_global(:), halo_global(:)
character(len=24) :: name
info = psb_success_
name = 'psb_cd_nest_compose'
if (.not. allocated(desc_grid%descs)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='nested descriptor not allocated')
return
end if
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
ctxt = desc_grid%descs(1,1)%get_context()
! 1. field offsets in the global numbering
allocate(field_offset(n_fields+1), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
field_offset(1) = 0
do i_field = 1, n_fields
field_offset(i_field+1) = field_offset(i_field) + desc_grid%descs(1,i_field)%get_global_rows()
end do
! 2. local owned global indices: U_k { offset_k + l2g(owned of field i_field) }
owned_count = 0
do i_field = 1, n_fields
owned_count = owned_count + desc_grid%descs(1,i_field)%get_local_rows()
end do
allocate(owned_global(owned_count), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
owned_count = 0
do i_field = 1, n_fields
n_owned = desc_grid%descs(1,i_field)%get_local_rows()
do i_loc = 1, n_owned
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'); return
end if
owned_count = owned_count + 1
owned_global(owned_count) = field_offset(i_field) + global_idx
end do
end do
! 3. allocate the global descriptor with the concatenated ownership
call psb_cdall(ctxt, desc_global, info, vl=owned_global)
if (info /= 0) then
call psb_errpush(psb_err_from_subroutine_, name, a_err='psb_cdall'); return
end if
! 4. global halo: U_k { offset_k + l2g(halo of field i_field) }
! field-i_field halo local indices are local_rows+1 .. local_cols
halo_count = 0
do i_field = 1, n_fields
halo_count = halo_count + (desc_grid%descs(1,i_field)%get_local_cols() &
& - desc_grid%descs(1,i_field)%get_local_rows())
end do
if (halo_count > 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
base/modules/tools/psb_d_nest_builder_mod.F90
+368
View File
@@ -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
base/modules/tools/psb_d_nest_tools_mod.F90
+137 -214
View File
@@ -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
info = psb_success_
name = 'psb_spasb_nest'
dupl_mode = psb_dupl_add_
if (present(dupl)) dupl_mode = 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_)
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