diff --git a/CMakeLists.txt b/CMakeLists.txt index c3b44d192..4e5ab79a9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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() diff --git a/base/modules/tools/psb_c_nest_builder_mod.F90 b/base/modules/tools/psb_c_nest_builder_mod.F90 index 27838cb02..31676a53e 100644 --- a/base/modules/tools/psb_c_nest_builder_mod.F90 +++ b/base/modules/tools/psb_c_nest_builder_mod.F90 @@ -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. diff --git a/base/modules/tools/psb_d_nest_builder_mod.F90 b/base/modules/tools/psb_d_nest_builder_mod.F90 index fcbc5ed5a..dc7e96dcf 100644 --- a/base/modules/tools/psb_d_nest_builder_mod.F90 +++ b/base/modules/tools/psb_d_nest_builder_mod.F90 @@ -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. diff --git a/base/modules/tools/psb_s_nest_builder_mod.F90 b/base/modules/tools/psb_s_nest_builder_mod.F90 index 59686b3f2..8cdbe093a 100644 --- a/base/modules/tools/psb_s_nest_builder_mod.F90 +++ b/base/modules/tools/psb_s_nest_builder_mod.F90 @@ -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. diff --git a/base/modules/tools/psb_z_nest_builder_mod.F90 b/base/modules/tools/psb_z_nest_builder_mod.F90 index c2f92f9df..26dcbf82f 100644 --- a/base/modules/tools/psb_z_nest_builder_mod.F90 +++ b/base/modules/tools/psb_z_nest_builder_mod.F90 @@ -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. diff --git a/linsolve/CMakeLists.txt b/linsolve/CMakeLists.txt index 039ba29dc..6092a2817 100644 --- a/linsolve/CMakeLists.txt +++ b/linsolve/CMakeLists.txt @@ -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 diff --git a/prec/psb_d_nestedprec.f90 b/prec/psb_d_nestedprec.f90 index 1978a17c7..0446f0706 100644 --- a/prec/psb_d_nestedprec.f90 +++ b/prec/psb_d_nestedprec.f90 @@ -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 diff --git a/test/nested/CMakeLists.txt b/test/nested/CMakeLists.txt index 463d908fc..1bd2650d5 100644 --- a/test/nested/CMakeLists.txt +++ b/test/nested/CMakeLists.txt @@ -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=") +# 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}") \ No newline at end of file diff --git a/test/nested/README.md b/test/nested/README.md index 5234a4191..ef4f83e8e 100644 --- a/test/nested/README.md +++ b/test/nested/README.md @@ -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. diff --git a/test/nested/psb_d_nest_cg_test.F90 b/test/nested/psb_d_nest_cg_test.F90 index 1148ccb81..52eeed0af 100644 --- a/test/nested/psb_d_nest_cg_test.F90 +++ b/test/nested/psb_d_nest_cg_test.F90 @@ -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 diff --git a/test/nested/psb_d_nest_glob_test.F90 b/test/nested/psb_d_nest_glob_test.F90 index f12b3f7d1..df79ffb84 100644 --- a/test/nested/psb_d_nest_glob_test.F90 +++ b/test/nested/psb_d_nest_glob_test.F90 @@ -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 diff --git a/test/nested/psb_d_nest_mult_prec_cg_test.F90 b/test/nested/psb_d_nest_mult_prec_cg_test.F90 index 0ee363b47..6a33e9ff7 100644 --- a/test/nested/psb_d_nest_mult_prec_cg_test.F90 +++ b/test/nested/psb_d_nest_mult_prec_cg_test.F90 @@ -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 diff --git a/test/nested/psb_d_nest_petsc_stokes_test.F90 b/test/nested/psb_d_nest_petsc_stokes_test.F90 index 2cd1b913b..2d4ba353c 100644 --- a/test/nested/psb_d_nest_petsc_stokes_test.F90 +++ b/test/nested/psb_d_nest_petsc_stokes_test.F90 @@ -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 \ No newline at end of file +end program psb_d_nest_petsc_stokes_test diff --git a/test/nested/psb_d_nest_prec_cg_test.F90 b/test/nested/psb_d_nest_prec_cg_test.F90 index 539bb44b9..10b5c2550 100644 --- a/test/nested/psb_d_nest_prec_cg_test.F90 +++ b/test/nested/psb_d_nest_prec_cg_test.F90 @@ -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 diff --git a/test/nested/psb_d_nest_rect_test.F90 b/test/nested/psb_d_nest_rect_test.F90 index 27e621cf0..5f975758e 100644 --- a/test/nested/psb_d_nest_rect_test.F90 +++ b/test/nested/psb_d_nest_rect_test.F90 @@ -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= 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