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

Fix nested matrix lifecycle and distributed preconditioning

Browse Source
This commit is contained in:
jalmerol
2026-09-21 14:52:05 +02:00
parent 157d5ea6ce
commit 0c93e0781b
16 changed files with 589 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+5
View File
@@ -1407,3 +1407,8 @@ message(STATUS "CMAKE_INSTALL_MODULDIR: ${CMAKE_INSTALL_MODULDIR} - ${PSB_CMAKE_
# Unit tests targeting each function, argument, and branch of code
# add_mpi_test(initialize_mpi 2 initialize_mpi)
option(PSB_BUILD_NESTED_TESTS "Build the nested matrix tests" OFF)
if(PSB_BUILD_NESTED_TESTS)
find_package(MPI REQUIRED COMPONENTS Fortran)
add_subdirectory(test/nested)
endif()
base/modules/tools/psb_c_nest_builder_mod.F90
+44 -9
View File
@@ -64,7 +64,7 @@
module psb_c_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_penv_mod, only : psb_ctxt_type, psb_info, psb_sum
use psb_desc_mod, only : psb_desc_type
use psb_c_mat_mod, only : psb_cspmat_type
use psb_c_base_mat_mod, only : psb_c_base_sparse_mat
@@ -126,6 +126,12 @@ contains
name = 'psb_c_nest_op_init'
call psb_info(context, my_rank, num_procs)
if (num_procs <= 0 .or. size(field_sizes) == 0 .or. any(field_sizes <= 0)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid context or field sizes'); return
end if
call op%free(info)
if (info /= psb_success_) return
n_fields = size(field_sizes)
op%context = context
op%n_fields = n_fields
@@ -171,7 +177,12 @@ contains
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
if (n_entries < 0 .or. n_entries > size(entry_rows) .or. &
& n_entries > size(entry_cols) .or. n_entries > size(entry_vals)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid triplet count'); 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)
@@ -200,11 +211,15 @@ contains
class(psb_c_base_sparse_mat), intent(in), optional :: mold
type(psb_c_nest_base_mat) :: nest_operator
integer(psb_ipk_) :: n_fields, i_field, j_field
integer(psb_ipk_) :: n_fields, i_field, j_field, block_present
character(len=24) :: name
info = psb_success_
name = 'psb_c_nest_op_asb'
if (op%assembled .or. op%n_fields <= 0 .or. .not. allocated(op%block_buffer)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='operator not in build state'); return
end if
n_fields = op%n_fields
! 1) assemble the per-field descriptors (with the union halo accumulated in ins)
@@ -224,7 +239,19 @@ contains
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
! Presence is global: every rank must build the same block structure,
! including empty local blocks on ranks owning no entries of a field.
block_present = merge(1, 0, op%block_buffer(i_field,j_field)%n_entries > 0)
call psb_sum(op%context, block_present)
if (block_present > 0) then
if (.not. allocated(op%block_buffer(i_field,j_field)%entry_rows)) then
allocate(op%block_buffer(i_field,j_field)%entry_rows(0), &
& op%block_buffer(i_field,j_field)%entry_cols(0), &
& op%block_buffer(i_field,j_field)%entry_vals(0), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
end if
call psb_c_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, &
@@ -249,6 +276,7 @@ contains
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)
if (info /= psb_success_) return
end do
end do
@@ -294,17 +322,24 @@ contains
end do
end do
deallocate(op%block_buffer, stat=local_info)
if (info == psb_success_) info = 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
! The adapter is a non-owning view; release the blocks explicitly. Do not
! gate cleanup on assembled: an earlier assembly may have failed midway.
call op%a_glob%free()
call op%block_storage%free(local_info)
if (info == psb_success_) info = local_info
call op%desc_glob%free(local_info)
if (info == psb_success_) info = local_info
call op%grid_desc%free(local_info)
if (info == psb_success_) info = local_info
if (allocated(op%field_desc)) then
do i_field = 1, size(op%field_desc)
call op%field_desc(i_field)%free(local_info)
if (info == psb_success_) info = local_info
end do
deallocate(op%field_desc, stat=local_info)
if (info == psb_success_) info = local_info
end if
op%n_fields = 0
op%assembled = .false.
base/modules/tools/psb_d_nest_builder_mod.F90
+44 -9
View File
@@ -64,7 +64,7 @@
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_penv_mod, only : psb_ctxt_type, psb_info, psb_sum
use psb_desc_mod, only : psb_desc_type
use psb_d_mat_mod, only : psb_dspmat_type
use psb_d_base_mat_mod, only : psb_d_base_sparse_mat
@@ -126,6 +126,12 @@ contains
name = 'psb_d_nest_op_init'
call psb_info(context, my_rank, num_procs)
if (num_procs <= 0 .or. size(field_sizes) == 0 .or. any(field_sizes <= 0)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid context or field sizes'); return
end if
call op%free(info)
if (info /= psb_success_) return
n_fields = size(field_sizes)
op%context = context
op%n_fields = n_fields
@@ -171,7 +177,12 @@ contains
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
if (n_entries < 0 .or. n_entries > size(entry_rows) .or. &
& n_entries > size(entry_cols) .or. n_entries > size(entry_vals)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid triplet count'); 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)
@@ -200,11 +211,15 @@ contains
class(psb_d_base_sparse_mat), intent(in), optional :: mold
type(psb_d_nest_base_mat) :: nest_operator
integer(psb_ipk_) :: n_fields, i_field, j_field
integer(psb_ipk_) :: n_fields, i_field, j_field, block_present
character(len=24) :: name
info = psb_success_
name = 'psb_d_nest_op_asb'
if (op%assembled .or. op%n_fields <= 0 .or. .not. allocated(op%block_buffer)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='operator not in build state'); return
end if
n_fields = op%n_fields
! 1) assemble the per-field descriptors (with the union halo accumulated in ins)
@@ -224,7 +239,19 @@ contains
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
! Presence is global: every rank must build the same block structure,
! including empty local blocks on ranks owning no entries of a field.
block_present = merge(1, 0, op%block_buffer(i_field,j_field)%n_entries > 0)
call psb_sum(op%context, block_present)
if (block_present > 0) then
if (.not. allocated(op%block_buffer(i_field,j_field)%entry_rows)) then
allocate(op%block_buffer(i_field,j_field)%entry_rows(0), &
& op%block_buffer(i_field,j_field)%entry_cols(0), &
& op%block_buffer(i_field,j_field)%entry_vals(0), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
end if
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, &
@@ -249,6 +276,7 @@ contains
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)
if (info /= psb_success_) return
end do
end do
@@ -294,17 +322,24 @@ contains
end do
end do
deallocate(op%block_buffer, stat=local_info)
if (info == psb_success_) info = 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
! The adapter is a non-owning view; release the blocks explicitly. Do not
! gate cleanup on assembled: an earlier assembly may have failed midway.
call op%a_glob%free()
call op%block_storage%free(local_info)
if (info == psb_success_) info = local_info
call op%desc_glob%free(local_info)
if (info == psb_success_) info = local_info
call op%grid_desc%free(local_info)
if (info == psb_success_) info = local_info
if (allocated(op%field_desc)) then
do i_field = 1, size(op%field_desc)
call op%field_desc(i_field)%free(local_info)
if (info == psb_success_) info = local_info
end do
deallocate(op%field_desc, stat=local_info)
if (info == psb_success_) info = local_info
end if
op%n_fields = 0
op%assembled = .false.
base/modules/tools/psb_s_nest_builder_mod.F90
+44 -9
View File
@@ -64,7 +64,7 @@
module psb_s_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_penv_mod, only : psb_ctxt_type, psb_info, psb_sum
use psb_desc_mod, only : psb_desc_type
use psb_s_mat_mod, only : psb_sspmat_type
use psb_s_base_mat_mod, only : psb_s_base_sparse_mat
@@ -126,6 +126,12 @@ contains
name = 'psb_s_nest_op_init'
call psb_info(context, my_rank, num_procs)
if (num_procs <= 0 .or. size(field_sizes) == 0 .or. any(field_sizes <= 0)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid context or field sizes'); return
end if
call op%free(info)
if (info /= psb_success_) return
n_fields = size(field_sizes)
op%context = context
op%n_fields = n_fields
@@ -171,7 +177,12 @@ contains
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
if (n_entries < 0 .or. n_entries > size(entry_rows) .or. &
& n_entries > size(entry_cols) .or. n_entries > size(entry_vals)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid triplet count'); 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)
@@ -200,11 +211,15 @@ contains
class(psb_s_base_sparse_mat), intent(in), optional :: mold
type(psb_s_nest_base_mat) :: nest_operator
integer(psb_ipk_) :: n_fields, i_field, j_field
integer(psb_ipk_) :: n_fields, i_field, j_field, block_present
character(len=24) :: name
info = psb_success_
name = 'psb_s_nest_op_asb'
if (op%assembled .or. op%n_fields <= 0 .or. .not. allocated(op%block_buffer)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='operator not in build state'); return
end if
n_fields = op%n_fields
! 1) assemble the per-field descriptors (with the union halo accumulated in ins)
@@ -224,7 +239,19 @@ contains
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
! Presence is global: every rank must build the same block structure,
! including empty local blocks on ranks owning no entries of a field.
block_present = merge(1, 0, op%block_buffer(i_field,j_field)%n_entries > 0)
call psb_sum(op%context, block_present)
if (block_present > 0) then
if (.not. allocated(op%block_buffer(i_field,j_field)%entry_rows)) then
allocate(op%block_buffer(i_field,j_field)%entry_rows(0), &
& op%block_buffer(i_field,j_field)%entry_cols(0), &
& op%block_buffer(i_field,j_field)%entry_vals(0), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
end if
call psb_s_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, &
@@ -249,6 +276,7 @@ contains
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)
if (info /= psb_success_) return
end do
end do
@@ -294,17 +322,24 @@ contains
end do
end do
deallocate(op%block_buffer, stat=local_info)
if (info == psb_success_) info = 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
! The adapter is a non-owning view; release the blocks explicitly. Do not
! gate cleanup on assembled: an earlier assembly may have failed midway.
call op%a_glob%free()
call op%block_storage%free(local_info)
if (info == psb_success_) info = local_info
call op%desc_glob%free(local_info)
if (info == psb_success_) info = local_info
call op%grid_desc%free(local_info)
if (info == psb_success_) info = local_info
if (allocated(op%field_desc)) then
do i_field = 1, size(op%field_desc)
call op%field_desc(i_field)%free(local_info)
if (info == psb_success_) info = local_info
end do
deallocate(op%field_desc, stat=local_info)
if (info == psb_success_) info = local_info
end if
op%n_fields = 0
op%assembled = .false.
base/modules/tools/psb_z_nest_builder_mod.F90
+44 -9
View File
@@ -64,7 +64,7 @@
module psb_z_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_penv_mod, only : psb_ctxt_type, psb_info, psb_sum
use psb_desc_mod, only : psb_desc_type
use psb_z_mat_mod, only : psb_zspmat_type
use psb_z_base_mat_mod, only : psb_z_base_sparse_mat
@@ -126,6 +126,12 @@ contains
name = 'psb_z_nest_op_init'
call psb_info(context, my_rank, num_procs)
if (num_procs <= 0 .or. size(field_sizes) == 0 .or. any(field_sizes <= 0)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid context or field sizes'); return
end if
call op%free(info)
if (info /= psb_success_) return
n_fields = size(field_sizes)
op%context = context
op%n_fields = n_fields
@@ -171,7 +177,12 @@ contains
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
if (n_entries < 0 .or. n_entries > size(entry_rows) .or. &
& n_entries > size(entry_cols) .or. n_entries > size(entry_vals)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='invalid triplet count'); 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)
@@ -200,11 +211,15 @@ contains
class(psb_z_base_sparse_mat), intent(in), optional :: mold
type(psb_z_nest_base_mat) :: nest_operator
integer(psb_ipk_) :: n_fields, i_field, j_field
integer(psb_ipk_) :: n_fields, i_field, j_field, block_present
character(len=24) :: name
info = psb_success_
name = 'psb_z_nest_op_asb'
if (op%assembled .or. op%n_fields <= 0 .or. .not. allocated(op%block_buffer)) then
info = psb_err_invalid_input_
call psb_errpush(info, name, a_err='operator not in build state'); return
end if
n_fields = op%n_fields
! 1) assemble the per-field descriptors (with the union halo accumulated in ins)
@@ -224,7 +239,19 @@ contains
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
! Presence is global: every rank must build the same block structure,
! including empty local blocks on ranks owning no entries of a field.
block_present = merge(1, 0, op%block_buffer(i_field,j_field)%n_entries > 0)
call psb_sum(op%context, block_present)
if (block_present > 0) then
if (.not. allocated(op%block_buffer(i_field,j_field)%entry_rows)) then
allocate(op%block_buffer(i_field,j_field)%entry_rows(0), &
& op%block_buffer(i_field,j_field)%entry_cols(0), &
& op%block_buffer(i_field,j_field)%entry_vals(0), stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_; call psb_errpush(info, name); return
end if
end if
call psb_z_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, &
@@ -249,6 +276,7 @@ contains
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)
if (info /= psb_success_) return
end do
end do
@@ -294,17 +322,24 @@ contains
end do
end do
deallocate(op%block_buffer, stat=local_info)
if (info == psb_success_) info = 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
! The adapter is a non-owning view; release the blocks explicitly. Do not
! gate cleanup on assembled: an earlier assembly may have failed midway.
call op%a_glob%free()
call op%block_storage%free(local_info)
if (info == psb_success_) info = local_info
call op%desc_glob%free(local_info)
if (info == psb_success_) info = local_info
call op%grid_desc%free(local_info)
if (info == psb_success_) info = local_info
if (allocated(op%field_desc)) then
do i_field = 1, size(op%field_desc)
call op%field_desc(i_field)%free(local_info)
if (info == psb_success_) info = local_info
end do
deallocate(op%field_desc, stat=local_info)
if (info == psb_success_) info = local_info
end if
op%n_fields = 0
op%assembled = .false.
linsolve/CMakeLists.txt
+4
View File
@@ -4,14 +4,17 @@ set(PSB_linsolve_source_files
impl/psb_crichardson.f90
impl/psb_zfcg.F90
impl/psb_zkrylov.f90
impl/psb_zminres.f90
impl/psb_srichardson.f90
impl/psb_ckrylov.f90
impl/psb_cminres.f90
impl/psb_crgmres.f90
impl/psb_cgcr.f90
impl/psb_sbicg.f90
impl/psb_ccgstabl.f90
impl/psb_drgmres.f90
impl/psb_skrylov.f90
impl/psb_sminres.f90
impl/psb_scg.F90
impl/psb_dcgstabl.f90
impl/psb_ccg.F90
@@ -26,6 +29,7 @@ set(PSB_linsolve_source_files
impl/psb_cfcg.F90
impl/psb_zrichardson.f90
impl/psb_dkrylov.f90
impl/psb_dminres.f90
impl/psb_ccgstab.f90
impl/psb_zcg.F90
impl/psb_dbicg.f90
prec/psb_d_nestedprec.f90
+56 -18
View File
@@ -887,6 +887,7 @@ contains
return
end if
ybuf(:) = dzero
if (beta /= dzero) ybuf(1:y%get_nrows()) = y%get_vect()
call psb_d_nested_apply(alpha, prec, xbuf, beta, ybuf, desc_data, info, trans)
if (info == psb_success_) call y%set(ybuf(1:y%get_nrows()))
deallocate(ybuf)
@@ -1269,6 +1270,7 @@ contains
bnorm = sqrt(max(dzero, psb_gedot(rhs, rhs, field_desc, info)))
if (info /= psb_success_) goto 100
if (bnorm == dzero) goto 100
r(:) = rhs(:) / bnorm
z(:) = dzero
call prec%blocks(field)%pc%apply(done, r, dzero, z, field_desc, info, trans='N', work=wrk)
@@ -1283,25 +1285,34 @@ contains
if (info /= psb_success_) exit
denom = psb_gedot(p, q, field_desc, info)
if (info /= psb_success_) exit
if (abs(denom) <= epsilon(done)) exit
if (abs(denom) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
alpha = rz / denom
sol(:) = sol(:) + alpha * p(:)
r(:) = r(:) - alpha * q(:)
rnorm = sqrt(max(dzero, psb_gedot(r, r, field_desc, info)))
if (info /= psb_success_) exit
if (rnorm <= kctx%tol * bnorm) exit
if (rnorm <= kctx%tol) exit
z(:) = dzero
call prec%blocks(field)%pc%apply(done, r, dzero, z, field_desc, info, trans='N', work=wrk)
if (info /= psb_success_) exit
rz_new = psb_gedot(r, z, field_desc, info)
if (info /= psb_success_) exit
if (abs(rz) <= epsilon(done)) exit
if (abs(rz) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
beta = rz_new / rz
p(:) = z(:) + beta * p(:)
rz = rz_new
end do
100 continue
sol(:) = bnorm * sol(:)
deallocate(r, z, p, q, wrk)
end subroutine psb_d_nested_inner_cg
@@ -1346,15 +1357,25 @@ contains
bnorm = sqrt(max(dzero, psb_gedot(rhs, rhs, field_desc, info)))
if (info /= psb_success_) goto 100
if (bnorm == dzero) goto 100
r(:) = rhs(:) / bnorm
r0(:) = r(:)
do k = 1, max(1, kctx%itmax)
rho = psb_gedot(r0, r, field_desc, info)
if (info /= psb_success_) exit
if (abs(rho) <= epsilon(done)) exit
if (abs(rho) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
if (k == 1) then
p(:) = r(:)
else
if (abs(omega) <= epsilon(done)) exit
if (abs(omega) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
beta = (rho / rho_old) * (alpha / omega)
p(:) = r(:) + beta * (p(:) - omega * v(:))
end if
@@ -1367,12 +1388,16 @@ contains
if (info /= psb_success_) exit
denom = psb_gedot(r0, v, field_desc, info)
if (info /= psb_success_) exit
if (abs(denom) <= epsilon(done)) exit
if (abs(denom) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
alpha = rho / denom
s(:) = r(:) - alpha * v(:)
rnorm = sqrt(max(dzero, psb_gedot(s, s, field_desc, info)))
if (info /= psb_success_) exit
if (rnorm <= kctx%tol * bnorm) then
if (rnorm <= kctx%tol) then
sol(:) = sol(:) + alpha * ph(:)
exit
end if
@@ -1385,18 +1410,23 @@ contains
if (info /= psb_success_) exit
denom = psb_gedot(t, t, field_desc, info)
if (info /= psb_success_) exit
if (abs(denom) <= epsilon(done)) exit
if (abs(denom) <= tiny(done)) then
info = psb_err_invalid_mat_state_
call psb_errpush(info, 'nested inner solve', a_err='Krylov breakdown')
exit
end if
omega = psb_gedot(t, s, field_desc, info) / denom
if (info /= psb_success_) exit
sol(:) = sol(:) + alpha * ph(:) + omega * sh(:)
r(:) = s(:) - omega * t(:)
rnorm = sqrt(max(dzero, psb_gedot(r, r, field_desc, info)))
if (info /= psb_success_) exit
if (rnorm <= kctx%tol * bnorm) exit
if (rnorm <= kctx%tol) exit
rho_old = rho
end do
100 continue
sol(:) = bnorm * sol(:)
deallocate(r, r0, p, v, s, t, ph, sh, wrk)
end subroutine psb_d_nested_inner_bicgstab
@@ -1456,7 +1486,7 @@ contains
block12 => psb_d_nest_get_block(prec%nest_op, 1, 2)
block21 => psb_d_nest_get_block(prec%nest_op, 2, 1)
if (associated(block12) .and. associated(block21)) then
block ! All ranks participate; absent local blocks contribute zero.
call psb_d_nest_apply_block(prec%nest_op, 1, 2, done, x2h, dzero, t1, info)
if (info /= psb_success_) goto 100
call psb_d_nested_field_solve(prec, 1, t1, w1, info)
@@ -1473,7 +1503,7 @@ contains
call psb_d_nest_apply_block(prec%nest_op, 2, 1, done, w1h, dzero, t2, info)
if (info /= psb_success_) goto 100
y2(:) = y2(:) - t2(:)
end if
end block
100 continue
deallocate(x2h, t1, w1, w1h, t2)
@@ -1752,7 +1782,9 @@ contains
if (info /= psb_success_) exit
res(:) = rhs(:) - sx(:)
if (prec%schur_tol > dzero) then
rnrm = sqrt(sum(res(1:n_owned) * res(1:n_owned)))
desc2 => psb_d_nest_get_field_desc(prec%nest_op, 2)
rnrm = sqrt(max(dzero, psb_gedot(res, res, desc2, info)))
if (info /= psb_success_) exit
if (rnrm <= prec%schur_tol) exit
end if
call psb_d_nested_field_solve(prec, 2, res, dz, info)
@@ -1948,7 +1980,7 @@ contains
block32 => psb_d_nest_get_block(prec%nest_op, 3, 2)
block33 => psb_d_nest_get_block(prec%nest_op, 3, 3)
if (associated(block13) .and. associated(block31)) then
block ! All ranks participate; absent local blocks contribute zero.
call psb_d_nest_apply_block(prec%nest_op, 1, 3, done, x3h, dzero, t1, info)
if (info /= psb_success_) goto 100
call psb_d_nested_field_solve(prec, 1, t1, w1, info)
@@ -1964,9 +1996,9 @@ contains
call psb_d_nest_apply_block(prec%nest_op, 3, 1, done, w1h, done, y3, info)
if (info /= psb_success_) goto 100
end if
end block
if (associated(block23) .and. associated(block32)) then
block ! All ranks participate; absent local blocks contribute zero.
call psb_d_nest_apply_block(prec%nest_op, 2, 3, done, x3h, dzero, t2, info)
if (info /= psb_success_) goto 100
call psb_d_nested_field_solve(prec, 2, t2, w2, info)
@@ -1982,7 +2014,7 @@ contains
call psb_d_nest_apply_block(prec%nest_op, 3, 2, done, w2h, done, y3, info)
if (info /= psb_success_) goto 100
end if
end block
if (associated(block33)) then
call psb_d_nest_apply_block(prec%nest_op, 3, 3, done, x3h, dzero, t3, info)
@@ -2041,8 +2073,13 @@ contains
r(:) = q(:)
rhs_norm = sqrt(max(dzero, psb_gedot(q, q, desc3, info)))
if (info /= psb_success_) goto 100
tol = prec%schur_tol
if (tol <= dzero) tol = 1.0e-6_psb_dpk_ * max(done, rhs_norm)
if (rhs_norm == dzero) goto 100
! Normalize to make the default relative tolerance independent of RHS scale.
q(:) = q(:) / rhs_norm
r(:) = q(:)
tol = 1.0e-6_psb_dpk_
if (prec%schur_tol > dzero) tol = prec%schur_tol / rhs_norm
if (done <= tol) goto 100
call psb_d_nested_field_solve(prec, 3, r, zc, info)
if (info /= psb_success_) goto 100
@@ -2086,6 +2123,7 @@ contains
end do
100 continue
sol(:) = rhs_norm * sol(:)
deallocate(q, r, zc, p, ap)
end subroutine psb_d_nested_pde_control_schur_solve
test/nested/CMakeLists.txt
+32 -39
View File
@@ -1,62 +1,55 @@
cmake_minimum_required(VERSION 3.10)
project(nested Fortran)
cmake_minimum_required(VERSION 3.15)
# Check for the installation path for psblas
if(NOT DEFINED PSBLAS_INSTALL_DIR)
message(FATAL_ERROR "Please specify the path to the psblas installation directory using -DPSBLAS_INSTALL_DIR=<path>")
# This directory can also be built independently against an installed PSBLAS.
if(CMAKE_SOURCE_DIR STREQUAL CMAKE_CURRENT_SOURCE_DIR)
project(PSBLASNestedTests Fortran)
find_package(psblas REQUIRED HINTS ${PSBLAS_INSTALL_DIR})
find_package(MPI REQUIRED COMPONENTS Fortran)
set(NEST_LIBS psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base MPI::MPI_Fortran)
if(PSBLAS_INSTALL_DIR)
include_directories("${PSBLAS_INSTALL_DIR}/include" "${PSBLAS_INSTALL_DIR}/modules")
endif()
else()
set(NEST_LIBS util linsolve prec ext base MPI::MPI_Fortran)
include_directories("${PROJECT_BINARY_DIR}/include" "${PROJECT_BINARY_DIR}/modules")
endif()
# Set the include and library directories based on the provided path
set(INSTALLDIR "${PSBLAS_INSTALL_DIR}")
set(INCDIR "${INSTALLDIR}/include")
set(MODDIR "${INSTALLDIR}/modules")
set(LIBDIR "${INSTALLDIR}/lib")
# Find the psblas package
find_package(psblas REQUIRED PATHS ${INSTALLDIR})
# Include directories for the Fortran compiler
include_directories(${INCDIR} ${MODDIR})
# Define executable directory
set(EXEDIR "${CMAKE_CURRENT_SOURCE_DIR}/runs")
file(MAKE_DIRECTORY ${EXEDIR})
# Executable directory
set(EXEDIR "${CMAKE_CURRENT_BINARY_DIR}/runs")
# Nested (block-structured / MATNEST) tests
set(SOURCES_D_NEST_GLOB_TEST psb_d_nest_glob_test.F90)
set(SOURCES_D_NEST_RECT_TEST psb_d_nest_rect_test.F90)
set(SOURCES_D_NEST_CG_TEST psb_d_nest_cg_test.F90)
set(SOURCES_D_NEST_GLOB_TEST psb_d_nest_glob_test.F90)
set(SOURCES_D_NEST_RECT_TEST psb_d_nest_rect_test.F90)
set(SOURCES_D_NEST_CG_TEST psb_d_nest_cg_test.F90)
set(SOURCES_D_NEST_PREC_CG_TEST psb_d_nest_prec_cg_test.F90)
set(SOURCES_D_NEST_MULT_PREC_CG_TEST psb_d_nest_mult_prec_cg_test.F90)
set(SOURCES_D_NEST_SCHUR_PREC_CG_TEST psb_d_nest_schur_prec_cg_test.F90)
set(SOURCES_D_NEST_PETSC_STOKES_TEST psb_d_nest_petsc_stokes_test.F90)
add_executable(psb_d_nest_glob_test ${SOURCES_D_NEST_GLOB_TEST})
target_link_libraries(psb_d_nest_glob_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_glob_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_glob_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_rect_test ${SOURCES_D_NEST_RECT_TEST})
target_link_libraries(psb_d_nest_rect_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_rect_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_rect_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_cg_test ${SOURCES_D_NEST_CG_TEST})
target_link_libraries(psb_d_nest_cg_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_cg_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_cg_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_prec_cg_test ${SOURCES_D_NEST_PREC_CG_TEST})
target_link_libraries(psb_d_nest_prec_cg_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_prec_cg_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_prec_cg_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_mult_prec_cg_test ${SOURCES_D_NEST_MULT_PREC_CG_TEST})
target_link_libraries(psb_d_nest_mult_prec_cg_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_mult_prec_cg_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_mult_prec_cg_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_schur_prec_cg_test ${SOURCES_D_NEST_SCHUR_PREC_CG_TEST})
target_link_libraries(psb_d_nest_schur_prec_cg_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
target_link_libraries(psb_d_nest_schur_prec_cg_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_schur_prec_cg_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
add_executable(psb_d_nest_petsc_stokes_test ${SOURCES_D_NEST_PETSC_STOKES_TEST})
target_link_libraries(psb_d_nest_petsc_stokes_test psblas::util psblas::linsolve psblas::prec psblas::ext psblas::base)
# Set output directory for executables
foreach(target psb_d_nest_glob_test psb_d_nest_rect_test psb_d_nest_cg_test psb_d_nest_prec_cg_test psb_d_nest_mult_prec_cg_test psb_d_nest_schur_prec_cg_test psb_d_nest_petsc_stokes_test)
set_target_properties(${target} PROPERTIES
RUNTIME_OUTPUT_DIRECTORY ${EXEDIR}
)
endforeach()
target_link_libraries(psb_d_nest_petsc_stokes_test PRIVATE ${NEST_LIBS})
set_target_properties(psb_d_nest_petsc_stokes_test PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${EXEDIR}")
test/nested/README.md
+18
View File
@@ -273,3 +273,21 @@ Library (under `base/modules/`):
* `tools/psb_d_nest_tools_mod.F90` — block tools (`psb_d_nest_rect_block`, ...)
* `tools/psb_d_nest_builder_mod.F90` — `psb_d_nest_matrix` frontend (init/ins/asb)
* `psb_d_nest_mod.f90` — umbrella module (`use psb_d_nest_mod`)
## Automated nested tests
Standalone CMake builds register the six synthetic nested tests on 1, 2, 3,
and 4 MPI ranks: global and rectangular matrix products, CG, additive and
multiplicative preconditioning, and Schur preconditioning.
```sh
cmake -S test/nested -B build-nested -DPSBLAS_INSTALL_DIR=/path/to/psblas
cmake --build build-nested -j2
ctest --test-dir build-nested --output-on-failure
```
Alternatively configure the main project with `-DPSB_BUILD_NESTED_TESTS=ON`.
`NEST_TEST_RANKS` controls MPI sizes. Set `PETSC_STOKES_DIR` at configure time
to register the external PETSc binary fixture tests. Failed numerical tolerance
checks and library errors cause nonzero exits; CTest timeouts also catch MPI
hangs.
test/nested/psb_d_nest_cg_test.F90
+36 -2
View File
@@ -75,7 +75,7 @@ program psb_d_nest_cg_test
integer(psb_ipk_) :: field1_local_rows, field2_local_rows
integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dprec_type) :: preconditioner
type(psb_d_vect_type) :: x_solution, rhs, x_exact
@@ -108,6 +108,7 @@ program psb_d_nest_cg_test
! 1) create the nested operator: 2 fields of global size field_size
!---------------------------------------------------------------
call nested_matrix%init(context, [field_size, field_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%init info=', info; goto 9999
end if
@@ -128,6 +129,7 @@ program psb_d_nest_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(1, 1, field1_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,2) = diag*I
@@ -139,6 +141,7 @@ program psb_d_nest_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(2, 2, field2_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (1,2) = C : rows field1, cols field2 ; C(r,r)=-1, C(r,r-1)=-1
@@ -158,6 +161,7 @@ program psb_d_nest_cg_test
end if
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,1) = C^T : rows field2, cols field1 ; C^T(s,s)=-1, C^T(s,s+1)=-1
@@ -177,12 +181,14 @@ program psb_d_nest_cg_test
end if
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
!---------------------------------------------------------------
! 3) assemble: nested_matrix%a_glob / nested_matrix%desc_glob are ready for Krylov
!---------------------------------------------------------------
call nested_matrix%asb(info)
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%asb info=', info; goto 9999
end if
@@ -191,16 +197,23 @@ program psb_d_nest_cg_test
! 4) consistent RHS: x_exact = 1, rhs = M * x_exact (via the nested operator)
!---------------------------------------------------------------
call psb_geall(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call x_exact%set(done) ! x_exact = 1 everywhere
call psb_geall(rhs, nested_matrix%desc_glob, info); call psb_geasb(rhs, nested_matrix%desc_glob, info)
call psb_geall(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_exact, dzero, rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (RHS) info=', info
goto 9999
end if
norm_x_exact = psb_genrm2(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
!---------------------------------------------------------------
! 5) solve with the standard PSBLAS CG under every stock preconditioner
@@ -211,17 +224,22 @@ program psb_d_nest_cg_test
iter_diag = -1
do i_prec = 1, n_precs
call preconditioner%init(context, trim(prec_names(i_prec)), info)
call check_info(info, 'preconditioner%init')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: prec%build (', trim(prec_names(i_prec)), ') info=', info
all_passed = .false.; exit
end if
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG,', trim(prec_names(i_prec)), ') info=', info
all_passed = .false.; exit
@@ -229,7 +247,9 @@ program psb_d_nest_cg_test
! solution error: || x_solution - x_exact || / || x_exact ||
call psb_geaxpby(-done, x_exact, done, x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info) / norm_x_exact
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,a6,a,i6,a,es12.4,a,es12.4)') ' prec=', prec_names(i_prec), &
@@ -241,7 +261,9 @@ program psb_d_nest_cg_test
if (trim(prec_names(i_prec)) == 'DIAG') iter_diag = n_iter
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
end do
!---------------------------------------------------------------
@@ -254,12 +276,24 @@ program psb_d_nest_cg_test
write(*,*) '[PASS] CG converges on the nested operator with NONE/DIAG/BJAC'
else
write(*,*) '[FAIL] preconditioned CG on the nested operator (tol ', solution_tol, ')'
call psb_abort(context)
end if
end if
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_cg_test
test/nested/psb_d_nest_glob_test.F90
+41 -3
View File
@@ -64,7 +64,7 @@ program psb_d_nest_glob_test
integer(psb_ipk_) :: entry_idx, field1_local_rows, field2_local_rows
integer(psb_lpk_) :: global_row, global_col, field_size
type(psb_d_nest_matrix) :: nested_matrix ! the nested operator (init/ins/asb)
type(psb_d_nest_matrix), target :: nested_matrix ! the nested operator (init/ins/asb)
type(psb_dspmat_type) :: monolithic_ref ! monolithic CSR oracle
type(psb_d_vect_type) :: x_vec, y_nested, y_monolithic
@@ -84,6 +84,7 @@ program psb_d_nest_glob_test
! 1) build the 2x2 nested operator through the utility
!---------------------------------------------------------------
call nested_matrix%init(context, [field_size, field_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%init info=', info; goto 9999
end if
@@ -119,6 +120,7 @@ program psb_d_nest_glob_test
end if
end do
call nested_matrix%ins(1, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! B^T = 0.5 I -> block (1,2): rows in field 1, columns in field 2
@@ -132,6 +134,7 @@ program psb_d_nest_glob_test
entry_vals(entry_idx) = 0.5_psb_dpk_
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! B = 0.3 I -> block (2,1): rows in field 2, columns in field 1
@@ -145,11 +148,13 @@ program psb_d_nest_glob_test
entry_vals(entry_idx) = 0.3_psb_dpk_
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! assemble with the blocks stored in HLL (psb_ext format): exercises the
! configurable block storage and the format-agnostic nested matvec
call nested_matrix%asb(info, mold=hll_mold)
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%asb info=', info; goto 9999
end if
@@ -160,53 +165,74 @@ program psb_d_nest_glob_test
!---------------------------------------------------------------
call psb_spall(monolithic_ref, nested_matrix%desc_glob, info, &
& nnz=5*nested_matrix%desc_glob%get_local_rows())
call check_info(info, 'psb_spall')
do i_local_row = 1, field1_local_rows ! field-1 rows
global_row = field1_rows(i_local_row)
insert_value(1) = 2.0_psb_dpk_
call psb_spins(1,[global_row],[global_row],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
if (global_row > 1) then
insert_value(1)=-1.0_psb_dpk_
call psb_spins(1,[global_row],[global_row-1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end if
if (global_row < field_size) then
insert_value(1)=-1.0_psb_dpk_
call psb_spins(1,[global_row],[global_row+1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end if
global_col = field_size + global_row
insert_value(1) = 0.5_psb_dpk_ ! B^T
call psb_spins(1,[global_row],[global_col],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end do
do i_local_row = 1, field2_local_rows ! field-2 rows
global_row = field2_rows(i_local_row)
global_col = global_row
insert_value(1) = 0.3_psb_dpk_ ! B
call psb_spins(1,[field_size+global_row],[global_col],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end do
call psb_spasb(monolithic_ref, nested_matrix%desc_glob, info, dupl=psb_dupl_add_)
call check_info(info, 'psb_spasb')
!---------------------------------------------------------------
! 4) compare the two matrix-vector products on a distinct-valued x (x[g] = g)
!---------------------------------------------------------------
call psb_geall(x_vec, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
do i_local_row = 1, nested_matrix%desc_glob%get_local_rows()
call nested_matrix%desc_glob%l2g(i_local_row, global_row, info)
call check_info(info, 'nested_matrix%desc_glob%l2g')
insert_value(1) = real(global_row, psb_dpk_)
call psb_geins(1, [global_row], insert_value, x_vec, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geins')
end do
call psb_geasb(x_vec, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_geall(y_nested, nested_matrix%desc_glob, info); call psb_geasb(y_nested, nested_matrix%desc_glob, info)
call psb_geall(y_monolithic, nested_matrix%desc_glob, info); call psb_geasb(y_monolithic, nested_matrix%desc_glob, info)
call psb_geall(y_nested, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(y_nested, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_geall(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_vec, dzero, y_nested, nested_matrix%desc_glob, info) ! via nested csmv
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (nested) info=', info
goto 9999
end if
call psb_spmm(done, monolithic_ref, x_vec, dzero, y_monolithic, nested_matrix%desc_glob, info) ! CSR oracle
call check_info(info, 'psb_spmm')
call psb_geaxpby(done, y_nested, -done, y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
mismatch_norm = psb_genrm2(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(field)=', field_size
@@ -215,12 +241,24 @@ program psb_d_nest_glob_test
write(*,*) '[PASS] nested global operator matches monolithic CSR'
else
write(*,*) '[FAIL] mismatch above tolerance ', tolerance
call psb_abort(context)
end if
end if
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_glob_test
test/nested/psb_d_nest_mult_prec_cg_test.F90
+43 -1
View File
@@ -64,7 +64,7 @@ program psb_d_nest_mult_prec_cg_test
integer(psb_ipk_) :: field1_local_rows, field2_local_rows
integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dprec_type) :: preconditioner
type(psb_d_vect_type) :: x_solution, rhs, x_exact
@@ -95,6 +95,7 @@ program psb_d_nest_mult_prec_cg_test
all_passed = .true.
call nested_matrix%init(context, [field_size, field_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%init info=', info
goto 9999
@@ -113,6 +114,7 @@ program psb_d_nest_mult_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(1, 1, field1_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,2) = diag*I
@@ -124,6 +126,7 @@ program psb_d_nest_mult_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(2, 2, field2_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (1,2) = C
@@ -143,6 +146,7 @@ program psb_d_nest_mult_prec_cg_test
end if
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,1) = C^T
@@ -162,37 +166,50 @@ program psb_d_nest_mult_prec_cg_test
end if
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
call nested_matrix%asb(info)
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%asb info=', info
goto 9999
end if
call psb_geall(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call x_exact%set(done)
call psb_geall(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_exact, dzero, rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (RHS) info=', info
goto 9999
end if
norm_x_exact = psb_genrm2(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(global)=', 2*field_size
! Baseline: no preconditioner.
call preconditioner%init(context, 'NONE', info)
call check_info(info, 'preconditioner%init')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter_none, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG,NONE) info=', info
all_passed = .false.
@@ -200,13 +217,19 @@ program psb_d_nest_mult_prec_cg_test
if (my_rank == 0) write(*,'(a,i6,a,es12.4)') &
& ' prec=NONE CG iterations=', n_iter_none, ' residual=', final_residual
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
do iprec = 1, size(composition_names)
call preconditioner%init(context, 'NEST', info)
call check_info(info, 'preconditioner%init')
call preconditioner%set('COMPOSITION', trim(composition_names(iprec)), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info)
call check_info(info, 'preconditioner%set')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested prec%build info=', info, &
& ' composition=', trim(composition_names(iprec))
@@ -215,10 +238,13 @@ program psb_d_nest_mult_prec_cg_test
end if
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov(trim(method_names(iprec)), nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(', trim(method_names(iprec)), ',NEST) info=', info, &
& ' composition=', trim(composition_names(iprec))
@@ -226,7 +252,9 @@ program psb_d_nest_mult_prec_cg_test
end if
call psb_geaxpby(-done, x_exact, done, x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info) / norm_x_exact
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,a32,a,a16,a,i6,a,es12.4,a,es12.4)') &
& ' prec=NEST/', trim(composition_names(iprec)), &
@@ -238,7 +266,9 @@ program psb_d_nest_mult_prec_cg_test
if ((n_iter >= max_iter) .or. (solution_error > solution_tol)) all_passed = .false.
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
end do
if (my_rank == 0) then
@@ -246,12 +276,24 @@ program psb_d_nest_mult_prec_cg_test
write(*,*) '[PASS] Krylov solvers converge with multiplicative NEST preconditioners'
else
write(*,*) '[FAIL] Krylov solvers with multiplicative NEST preconditioners'
call psb_abort(context)
end if
end if
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_mult_prec_cg_test
test/nested/psb_d_nest_petsc_stokes_test.F90
+21 -5
View File
@@ -35,7 +35,7 @@ program psb_d_nest_petsc_stokes_test
end type petsc_vector
type(psb_ctxt_type) :: context
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dprec_type) :: preconditioner
type(psb_d_vect_type) :: rhs, x_solution, residual
@@ -153,11 +153,17 @@ program psb_d_nest_petsc_stokes_test
call check_info(info, 'prec%init')
if (psb_toupper(trim(ptype)) == 'NEST') then
call preconditioner%set('COMPOSITION', trim(composition), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SCHUR_SOLVE', trim(schur_solve), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SCHUR_MAXIT', schur_maxit, info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', trim(block_solve), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_SOLVE', trim(sub_solve), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_FILLIN', sub_fillin, info)
call check_info(info, 'preconditioner%set')
end if
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'prec%build')
@@ -171,12 +177,17 @@ program psb_d_nest_petsc_stokes_test
t_solve = psb_wtime() - t0
call psb_geall(residual, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(residual, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_geaxpby(done, rhs, dzero, residual, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
call psb_spmm(-done, nested_matrix%a_glob, x_solution, done, residual, nested_matrix%desc_glob, info)
call check_info(info, 'residual spmm')
residual_norm = psb_genrm2(residual, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
rhs_norm = psb_genrm2(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,i0,a,i0,a,i0,a,i0)') ' block sizes: n_u=', n_u, ' n_p=', n_p, ' total=', n_total, ' np=', num_procs
@@ -194,10 +205,15 @@ program psb_d_nest_petsc_stokes_test
9999 continue
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
call psb_gefree(residual, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call psb_gefree(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
call psb_exit(context)
contains
@@ -291,7 +307,7 @@ contains
stop 2
end if
read(unit) classid
if (classid /= petsc_mat_classid) stop 'Bad PETSc matrix classid'
if (classid /= petsc_mat_classid) error stop 'Bad PETSc matrix classid'
read(unit) nrow32
read(unit) ncol32
read(unit) nnz32
@@ -318,7 +334,7 @@ contains
mat%val(k) = real(vals64(k), psb_dpk_)
end do
end do
if (k /= mat%nnz) stop 'Bad PETSc matrix row lengths'
if (k /= mat%nnz) error stop 'Bad PETSc matrix row lengths'
end subroutine read_petsc_matrix
subroutine read_petsc_vector(filename, vec)
@@ -335,7 +351,7 @@ contains
stop 2
end if
read(unit) classid
if (classid /= petsc_vec_classid) stop 'Bad PETSc vector classid'
if (classid /= petsc_vec_classid) error stop 'Bad PETSc vector classid'
read(unit) n32
vec%n = int(n32, psb_lpk_)
allocate(vals64(vec%n), vec%val(vec%n))
@@ -443,4 +459,4 @@ contains
if (allocated(in_vals)) deallocate(in_vals)
end subroutine clear_triplets
end program psb_d_nest_petsc_stokes_test
end program psb_d_nest_petsc_stokes_test
test/nested/psb_d_nest_prec_cg_test.F90
+43 -1
View File
@@ -62,7 +62,7 @@ program psb_d_nest_prec_cg_test
integer(psb_ipk_) :: field1_local_rows, field2_local_rows
integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dprec_type) :: preconditioner
type(psb_d_vect_type) :: x_solution, rhs, x_exact
@@ -89,6 +89,7 @@ program psb_d_nest_prec_cg_test
all_passed = .true.
call nested_matrix%init(context, [field_size, field_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%init info=', info
goto 9999
@@ -107,6 +108,7 @@ program psb_d_nest_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(1, 1, field1_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,2) = diag*I
@@ -118,6 +120,7 @@ program psb_d_nest_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(2, 2, field2_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (1,2) = C
@@ -137,6 +140,7 @@ program psb_d_nest_prec_cg_test
end if
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,1) = C^T
@@ -156,37 +160,50 @@ program psb_d_nest_prec_cg_test
end if
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
call nested_matrix%asb(info)
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%asb info=', info
goto 9999
end if
call psb_geall(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call x_exact%set(done)
call psb_geall(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_exact, dzero, rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (RHS) info=', info
goto 9999
end if
norm_x_exact = psb_genrm2(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(global)=', 2*field_size
! Baseline: no preconditioner.
call preconditioner%init(context, 'NONE', info)
call check_info(info, 'preconditioner%init')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter_none, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG,NONE) info=', info
all_passed = .false.
@@ -194,13 +211,19 @@ program psb_d_nest_prec_cg_test
if (my_rank == 0) write(*,'(a,i6,a,es12.4)') &
& ' prec=NONE CG iterations=', n_iter_none, ' residual=', final_residual
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
! Specialized nested preconditioner: additive block-diagonal field split.
call preconditioner%init(context, 'NEST', info)
call check_info(info, 'preconditioner%init')
call preconditioner%set('COMPOSITION', 'ADDITIVE', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info)
call check_info(info, 'preconditioner%set')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested prec%build info=', info
all_passed = .false.
@@ -208,17 +231,22 @@ program psb_d_nest_prec_cg_test
end if
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG,NEST) info=', info
all_passed = .false.
end if
call psb_geaxpby(-done, x_exact, done, x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info) / norm_x_exact
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,i6,a,es12.4,a,es12.4)') &
& ' prec=NEST CG iterations=', n_iter, ' residual=', final_residual, &
@@ -233,14 +261,28 @@ program psb_d_nest_prec_cg_test
write(*,*) '[PASS] CG converges with the specialized NEST preconditioner'
else
write(*,*) '[FAIL] CG with specialized NEST preconditioner'
call psb_abort(context)
end if
end if
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_prec_cg_test
test/nested/psb_d_nest_rect_test.F90
+42 -3
View File
@@ -59,7 +59,7 @@ program psb_d_nest_rect_test
integer(psb_ipk_) :: entry_idx, v_local_rows, q_local_rows
integer(psb_lpk_) :: v_global_row, q_global_row, q_col, v_size, q_size
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dspmat_type) :: monolithic_ref
type(psb_d_vect_type) :: x_vec, y_nested, y_monolithic
real(psb_dpk_) :: insert_value(1)
@@ -80,6 +80,7 @@ program psb_d_nest_rect_test
! 1) build the 2x2 nested operator (fields V, Q)
!---------------------------------------------------------------
call nested_matrix%init(context, [v_size, q_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%init info=', info; goto 9999
end if
@@ -114,6 +115,7 @@ program psb_d_nest_rect_test
end if
end do
call nested_matrix%ins(1, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! B^T rectangular -> block (1,2), V x Q : row r -> col mod(r-1,nQ)+1, val 0.5
@@ -127,6 +129,7 @@ program psb_d_nest_rect_test
entry_vals(entry_idx) = 0.5_psb_dpk_
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! B rectangular -> block (2,1), Q x V : row q -> cols q and q+nQ, val 0.3
@@ -144,11 +147,13 @@ program psb_d_nest_rect_test
entry_vals(entry_idx) = 0.3_psb_dpk_
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! assemble with the blocks stored in CSC instead of the CSR default:
! exercises the configurable block storage on a base format
call nested_matrix%asb(info, type='CSC')
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL: nested_matrix%asb info=', info; goto 9999
end if
@@ -158,51 +163,73 @@ program psb_d_nest_rect_test
!---------------------------------------------------------------
call psb_spall(monolithic_ref, nested_matrix%desc_glob, info, &
& nnz=6*nested_matrix%desc_glob%get_local_rows())
call check_info(info, 'psb_spall')
do i_local_row = 1, v_local_rows ! V rows
v_global_row = v_rows(i_local_row)
insert_value(1)=2.0_psb_dpk_
call psb_spins(1,[v_global_row],[v_global_row],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
if (v_global_row>1) then
insert_value(1)=-1.0_psb_dpk_
call psb_spins(1,[v_global_row],[v_global_row-1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end if
if (v_global_row<v_size) then
insert_value(1)=-1.0_psb_dpk_
call psb_spins(1,[v_global_row],[v_global_row+1_psb_lpk_],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end if
q_col = v_size + (mod(v_global_row-1_psb_lpk_, q_size) + 1)
insert_value(1)=0.5_psb_dpk_ ! B^T
call psb_spins(1,[v_global_row],[q_col],insert_value,monolithic_ref,nested_matrix%desc_glob,info)
call check_info(info, 'psb_spins')
end do
do i_local_row = 1, q_local_rows ! Q rows
q_global_row = q_rows(i_local_row)
insert_value(1)=0.3_psb_dpk_
call psb_spins(1,[v_size+q_global_row],[q_global_row], insert_value,monolithic_ref,nested_matrix%desc_glob,info) ! col q
call check_info(info, 'psb_spins')
call psb_spins(1,[v_size+q_global_row],[q_global_row+q_size],insert_value,monolithic_ref,nested_matrix%desc_glob,info) ! col q+nQ
call check_info(info, 'psb_spins')
end do
call psb_spasb(monolithic_ref, nested_matrix%desc_glob, info, dupl=psb_dupl_add_)
call check_info(info, 'psb_spasb')
!---------------------------------------------------------------
! 4) compare the two matrix-vector products on x[g] = g
!---------------------------------------------------------------
call psb_geall(x_vec, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
do i_local_row = 1, nested_matrix%desc_glob%get_local_rows()
call nested_matrix%desc_glob%l2g(i_local_row, v_global_row, info)
call check_info(info, 'nested_matrix%desc_glob%l2g')
insert_value(1) = real(v_global_row, psb_dpk_)
call psb_geins(1, [v_global_row], insert_value, x_vec, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geins')
end do
call psb_geasb(x_vec, nested_matrix%desc_glob, info)
call psb_geall(y_nested, nested_matrix%desc_glob, info); call psb_geasb(y_nested, nested_matrix%desc_glob, info)
call psb_geall(y_monolithic, nested_matrix%desc_glob, info); call psb_geasb(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_geall(y_nested, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(y_nested, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_geall(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_vec, dzero, y_nested, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank==0) write(*,*) 'FAIL spmm nested info=', info; goto 9999
end if
call psb_spmm(done, monolithic_ref, x_vec, dzero, y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
call psb_geaxpby(done, y_nested, -done, y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
mismatch_norm = psb_genrm2(y_monolithic, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,i0,a,i0,a,i0)') ' np=', num_procs, ' |V|=', v_size, ' |Q|=', q_size
@@ -211,12 +238,24 @@ program psb_d_nest_rect_test
write(*,*) '[PASS] rectangular nested operator matches monolithic CSR'
else
write(*,*) '[FAIL] mismatch above tolerance ', tolerance
call psb_abort(context)
end if
end if
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_rect_test
test/nested/psb_d_nest_schur_prec_cg_test.F90
+72 -2
View File
@@ -69,7 +69,7 @@ program psb_d_nest_schur_prec_cg_test
integer(psb_ipk_) :: field1_local_rows, field2_local_rows
integer(psb_lpk_) :: field1_global_row, field2_global_row, field_size
type(psb_d_nest_matrix) :: nested_matrix
type(psb_d_nest_matrix), target :: nested_matrix
type(psb_dprec_type) :: preconditioner
type(psb_d_vect_type) :: x_solution, rhs, x_exact
@@ -100,6 +100,7 @@ program psb_d_nest_schur_prec_cg_test
all_passed = .true.
call nested_matrix%init(context, [field_size, field_size], info)
call check_info(info, 'nested_matrix%init')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%init info=', info
goto 9999
@@ -118,6 +119,7 @@ program psb_d_nest_schur_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(1, 1, field1_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,2) = diag*I
@@ -129,6 +131,7 @@ program psb_d_nest_schur_prec_cg_test
entry_vals(i_local_row) = diag_value
end do
call nested_matrix%ins(2, 2, field2_local_rows, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (1,2) = C
@@ -148,6 +151,7 @@ program psb_d_nest_schur_prec_cg_test
end if
end do
call nested_matrix%ins(1, 2, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
! block (2,1) = C^T
@@ -167,37 +171,50 @@ program psb_d_nest_schur_prec_cg_test
end if
end do
call nested_matrix%ins(2, 1, entry_idx, entry_rows, entry_cols, entry_vals, info)
call check_info(info, 'nested_matrix%ins')
deallocate(entry_rows, entry_cols, entry_vals)
call nested_matrix%asb(info)
call check_info(info, 'nested_matrix%asb')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested_matrix%asb info=', info
goto 9999
end if
call psb_geall(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call x_exact%set(done)
call psb_geall(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_spmm(done, nested_matrix%a_glob, x_exact, dzero, rhs, nested_matrix%desc_glob, info)
call check_info(info, 'psb_spmm')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_spmm (RHS) info=', info
goto 9999
end if
norm_x_exact = psb_genrm2(x_exact, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) write(*,'(a,i0,a,i0)') ' np=', num_procs, ' N(global)=', 2*field_size
! Baseline: no preconditioner.
call preconditioner%init(context, 'NONE', info)
call check_info(info, 'preconditioner%init')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('CG', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter_none, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(CG,NONE) info=', info
all_passed = .false.
@@ -205,27 +222,44 @@ program psb_d_nest_schur_prec_cg_test
if (my_rank == 0) write(*,'(a,i6,a,es12.4)') &
& ' prec=NONE CG iterations=', n_iter_none, ' residual=', final_residual
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
do iprec = 1, size(composition_names)
call preconditioner%init(context, 'NEST', info)
call check_info(info, 'preconditioner%init')
call preconditioner%set('COMPOSITION', trim(composition_names(iprec)), info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SCHUR_SOLVE', trim(schur_solve_names(iprec)), info)
call check_info(info, 'preconditioner%set')
if (trim(schur_solve_names(iprec)) == 'MATRIX_FREE') then
call preconditioner%set('SCHUR_MAXIT', 8, info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'BJAC', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'BJAC', info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'BJAC', info, idx=2)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_SOLVE', 'ILU', info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_SOLVE', 'ILU', info, idx=2)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_FILLIN', 0, info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SUB_FILLIN', 0, info, idx=2)
call check_info(info, 'preconditioner%set')
else
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info, idx=2)
call check_info(info, 'preconditioner%set')
end if
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested prec%build info=', info, &
& ' composition=', trim(composition_names(iprec))
@@ -234,10 +268,13 @@ program psb_d_nest_schur_prec_cg_test
end if
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call psb_krylov('BICGSTAB', nested_matrix%a_glob, preconditioner, rhs, x_solution, stop_tol, &
& nested_matrix%desc_glob, info, &
& itmax=max_iter, iter=n_iter, err=final_residual, itrace=trace_level, istop=stop_criterion)
call check_info(info, 'psb_krylov')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: psb_krylov(BICGSTAB,NEST) info=', info, &
& ' composition=', trim(composition_names(iprec))
@@ -245,7 +282,9 @@ program psb_d_nest_schur_prec_cg_test
end if
call psb_geaxpby(-done, x_exact, done, x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geaxpby')
solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info) / norm_x_exact
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,a32,a,a16,a,a16,a,i6,a,es12.4,a,es12.4)') &
& ' prec=NEST/', trim(composition_names(iprec)), &
@@ -258,21 +297,34 @@ program psb_d_nest_schur_prec_cg_test
if ((n_iter >= max_iter) .or. (solution_error > solution_tol)) all_passed = .false.
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
end do
! Exercise independent inner Krylov contexts inside the field solves.
call preconditioner%init(context, 'NEST', info)
call check_info(info, 'preconditioner%init')
call preconditioner%set('COMPOSITION', 'SCHUR_FULL', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('SCHUR_SOLVE', 'A22', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('BLOCK_SOLVE', 'DIAG', info)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_SOLVE', 'CG', info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_MAXIT', 10, info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_TOL', 1.0d-10, info, idx=1)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_SOLVE', 'CG', info, idx=2)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_MAXIT', 12, info, idx=2)
call check_info(info, 'preconditioner%set')
call preconditioner%set('INNER_TOL', 1.0d-11, info, idx=2)
call check_info(info, 'preconditioner%set')
call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%build')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: nested inner Krylov prec%build info=', info
all_passed = .false.
@@ -280,14 +332,18 @@ program psb_d_nest_schur_prec_cg_test
end if
call psb_geall(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geall')
call psb_geasb(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_geasb')
call preconditioner%apply(rhs, x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'preconditioner%apply')
if (info /= psb_success_) then
if (my_rank == 0) write(*,*) 'FAIL: NEST inner Krylov apply info=', info
all_passed = .false.
end if
solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'norm')
if (my_rank == 0) then
write(*,'(a,a32,a,a16,a,a16,a,es12.4)') &
& ' prec=NEST/', 'SCHUR_FULL', &
@@ -296,11 +352,13 @@ program psb_d_nest_schur_prec_cg_test
& ' ||P rhs||=', solution_error
end if
if ((info /= psb_success_) .or. (solution_error /= solution_error) .or. &
if ((info /= psb_success_) .or. &
& (solution_error <= dzero)) all_passed = .false.
call psb_gefree(x_solution, nested_matrix%desc_glob, info)
call check_info(info, 'psb_gefree')
call preconditioner%free(info)
call check_info(info, 'preconditioner%free')
! Transposed nested preconditioner application is intentionally unsupported.
! The implementation returns psb_err_transpose_not_n_unsupported_ for trans /= 'N'.
@@ -314,12 +372,24 @@ program psb_d_nest_schur_prec_cg_test
write(*,*) ' including per-field inner Krylov solves'
else
write(*,*) '[FAIL] Krylov solvers with Schur-style NEST preconditioners'
call psb_abort(context)
end if
end if
call nested_matrix%free(info)
call check_info(info, 'nested_matrix%free')
9999 continue
call psb_exit(context)
contains
subroutine check_info(status, label)
integer(psb_ipk_), intent(in) :: status
character(len=*), intent(in) :: label
if (status /= psb_success_) then
write(*,*) '[FAIL] ', trim(label), ' rank=', my_rank, ' info=', status
call psb_abort(context)
end if
end subroutine check_info
end program psb_d_nest_schur_prec_cg_test