added SCHUR_SOLVE=SELFP / SELF option and stokes problem test

update for self/selfp

test/nested/CMakeLists.txt

@@ -29,6 +29,7 @@ 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)
@@ -48,8 +49,13 @@ target_link_libraries(psb_d_nest_mult_prec_cg_test psblas::util psblas::linsolve
 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)
 
+
+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)
+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}
     )

test/nested/Makefile

@@ -6,7 +6,11 @@ include $(INCDIR)/Make.inc.psblas
 # Libraries used
 LIBDIR=$(INSTALLDIR)/lib
 PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_linsolve -lpsb_prec -lpsb_ext -lpsb_base
+ifeq ($(strip $(BLAS)),)
+LDLIBS= $(PSBLDLIBS) -lblas
+else
 LDLIBS= $(PSBLDLIBS)
+endif
 #
 # Compilers and such
 #
@@ -17,7 +21,7 @@ FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG).
 EXEDIR=./runs
 
 all: runsd psb_d_nest_glob_test psb_d_nest_rect_test psb_d_nest_cg_test psb_d_nest_prec_cg_test \
-	psb_d_nest_mult_prec_cg_test psb_d_nest_schur_prec_cg_test
+	psb_d_nest_mult_prec_cg_test psb_d_nest_schur_prec_cg_test psb_d_nest_petsc_stokes_test
 
 runsd:
 	(if test ! -d runs ; then mkdir runs; fi)
@@ -46,12 +50,19 @@ psb_d_nest_schur_prec_cg_test: psb_d_nest_schur_prec_cg_test.o
 	$(FLINK) psb_d_nest_schur_prec_cg_test.o -o psb_d_nest_schur_prec_cg_test $(PSBLAS_LIB) $(LDLIBS)
 	/bin/mv psb_d_nest_schur_prec_cg_test $(EXEDIR)
 
+
+psb_d_nest_petsc_stokes_test: psb_d_nest_petsc_stokes_test.o
+	$(FLINK) psb_d_nest_petsc_stokes_test.o -o psb_d_nest_petsc_stokes_test $(PSBLAS_LIB) $(LDLIBS)
+	/bin/mv psb_d_nest_petsc_stokes_test $(EXEDIR)
+
+
 clean:
 	/bin/rm -f psb_d_nest_glob_test.o psb_d_nest_rect_test.o psb_d_nest_cg_test.o \
-	psb_d_nest_prec_cg_test.o psb_d_nest_mult_prec_cg_test.o psb_d_nest_schur_prec_cg_test.o *$(.mod) \
+	psb_d_nest_prec_cg_test.o psb_d_nest_mult_prec_cg_test.o psb_d_nest_schur_prec_cg_test.o \
+	psb_d_nest_petsc_stokes_test.o *$(.mod) \
 	$(EXEDIR)/psb_d_nest_glob_test $(EXEDIR)/psb_d_nest_rect_test $(EXEDIR)/psb_d_nest_cg_test \
 	$(EXEDIR)/psb_d_nest_prec_cg_test $(EXEDIR)/psb_d_nest_mult_prec_cg_test \
-	$(EXEDIR)/psb_d_nest_schur_prec_cg_test
+	$(EXEDIR)/psb_d_nest_schur_prec_cg_test $(EXEDIR)/psb_d_nest_petsc_stokes_test
 verycleanlib:
 	(cd ../..; make veryclean)
 lib:

test/nested/Nested_Preconditioner_Report.md

@@ -527,9 +527,11 @@ y
 
 No Schur matrix is ever stored. The code only needs block products and field solves.
 
-### Richardson Schur Solve
+### Matrix-Free Schur Solves
 
-When `SCHUR_SOLVE = MATRIX_FREE`, the preconditioner needs an approximation to `S^{-1}` even though it only has a routine for `S x`. The implementation therefore uses a small preconditioned Richardson iteration:
+When `SCHUR_SOLVE = MATRIX_FREE`, the preconditioner needs an approximation to a Schur inverse even though it only has routines that apply Schur products. The implementation uses two matrix-free strategies, depending on the selected composition.
+
+For the two-field Schur compositions (`SCHUR_LOWER`, `SCHUR_UPPER`, and `SCHUR_FULL`), the matrix-free Schur solve uses a small preconditioned Richardson iteration:
 
 ```text
 z_{k+1} = z_k + M^{-1} (b - S z_k)
@@ -543,7 +545,15 @@ Each Richardson step uses:
 2. a Schur residual `b - S z_k`
 3. the field-2 preconditioner or field-2 inner solve to correct `z_k`
 
-This gives a cheap inner solve for the Schur complement without assembling or storing `S`.
+For the three-field PDE-control composition (`SCHUR_PDE_CONTROL`), the Schur field is field 3 and the KKT Schur complement is negative for the standard distributed-control block signs. The implementation therefore solves the positive matrix-free equation
+
+```text
+T x_3 = -r_3,       T = -S
+```
+
+with a small conjugate-gradient iteration. The product `T x_3` is applied without assembling `T`; each CG iteration calls the PDE-control Schur action and uses `psb_gedot` with the field-3 descriptor for global reductions. The field-3 block preconditioner is still used as an optional left preconditioner, but when `A_33 = 0` it is only an identity fallback.
+
+This gives a stronger Schur inverse approximation for KKT optimal-control systems where the diagonal field-3 block is zero and a Richardson update based on `A_33` is ineffective.
 
 The implementation does not assemble an exact sparse Schur complement. Instead, it applies Schur actions using:
 
@@ -556,13 +566,23 @@ Supported Schur compositions are:
 - `SCHUR_LOWER`
 - `SCHUR_UPPER`
 - `SCHUR_FULL`
+- `SCHUR_PDE_CONTROL` (aliases: `PDE_CONTROL_SCHUR`, `SCHUR_2PLUS1`, `SCHUR3`, `SCHUR_3FIELD`, `SCHUR_THREE_FIELD`)
 
 The Schur solve mode is controlled by `SCHUR_SOLVE`:
 
 - `A22`: use the field-2 block preconditioner as the Schur approximation
-- `MATRIX_FREE`: use a small Richardson iteration around the matrix-free Schur action
+- `MATRIX_FREE`: use a matrix-free Schur solve; two-field Schur compositions use Richardson, while `SCHUR_PDE_CONTROL` uses CG on the positive field-3 Schur action
+- `SELFP` or `SELF`: build a PETSc-style lumped Schur approximation from `A_21`, the diagonal of `A_11`, and `A_12`, then build a normal PSBLAS `BJAC` Schur preconditioner on that approximation
 
-The `MATRIX_FREE` mode is controlled by `SCHUR_MAXIT` and `SCHUR_TOL`.
+For two-field Schur compositions, `SELFP` approximates
+
+```text
+S_p ~= A_22 - A_21 diag(A_11)^{-1} A_12
+```
+
+The assembled `S_p` is not the exact Schur complement. It is a sparse local approximation intended for use as a preconditioner, analogous in spirit to PETSc field-split `selfp` with a lumped inverse for `A_11`. After assembly, the nested preconditioner attaches the existing block-preconditioner machinery to `S_p`: currently this is a `BJAC` block preconditioner, and routed sub-options such as `SUB_SOLVE=ILU` and `SUB_FILLIN` are replayed onto it. During Schur application, `SCHUR_SOLVE=SELFP` applies this built Schur preconditioner directly instead of running Richardson on the matrix-free Schur action.
+
+The `MATRIX_FREE` mode is controlled by `SCHUR_MAXIT` and `SCHUR_TOL`. The `SELFP` mode is controlled primarily by the sub-preconditioner options routed to the Schur approximation, especially `SUB_SOLVE` and `SUB_FILLIN`.
 
 ## Implementation Structure
 
@@ -795,6 +815,39 @@ S x_2 = A_22 x_2 - A_21 B_1 A_12 x_2
 
 where `B_1` is the field-1 block solve or inner field solve. `psb_d_nested_schur_solve` either applies the field-2 block solve directly (`A22`) or runs a short preconditioned Richardson iteration using the matrix-free Schur action (`MATRIX_FREE`).
 
+## PDE-Control Schur Apply Path
+
+`psb_d_nested_apply_pde_control_schur` is defined for exactly three fields. It targets KKT systems with the block structure
+
+```text
+[ A_11   0    A_13 ]
+[  0    A_22  A_23 ]
+[ A_31  A_32  A_33 ]
+```
+
+Fields 1 and 2 are grouped into the leading diagonal block `D = diag(A_11,A_22)`, and field 3 is the Schur field. The preconditioner applies the usual block factorization steps:
+
+1. solve the leading fields with `B_1` and `B_2`
+2. update the field-3 residual with `A_31 x_1 + A_32 x_2`
+3. approximately solve the field-3 Schur equation
+4. correct fields 1 and 2 with the upper coupling through `A_13` and `A_23`
+
+The matrix-free positive PDE-control Schur action is
+
+```text
+T x_3 = A_31 B_1 A_13 x_3 + A_32 B_2 A_23 x_3 - A_33 x_3
+```
+
+For a standard distributed-control KKT matrix such as
+
+```text
+[ M      0      K  ]
+[ 0    alpha M -M  ]
+[ K     -M      0  ]
+```
+
+this operator is `T = -S`, where `S` is the field-3 KKT Schur complement. `psb_d_nested_pde_control_schur_solve` solves `T x_3 = -r_3` with a matrix-free CG iteration controlled by `SCHUR_MAXIT` and `SCHUR_TOL`. This is important when `A_33 = 0`: the old field-3 diagonal-block approximation cannot precondition the Schur complement, but the matrix-free CG path can still use the actual Schur action.
+
 ## Public Options
 
 The nested preconditioner exposes integer option keys in `psb_d_nestedprec`:
@@ -816,9 +869,9 @@ The user-facing string options routed by `psb_d_prec_type_impl.f90` are:
 
 | Option | Type | Values | Scope |
 | --- | --- | --- | --- |
-| `COMPOSITION`, `NEST_COMPOSITION` | character | `ADDITIVE`, `DIAGONAL`, `MULTIPLICATIVE`, `SYMMETRIC_MULTIPLICATIVE`, `SCHUR_LOWER`, `SCHUR_UPPER`, `SCHUR_FULL` | global |
+| `COMPOSITION`, `NEST_COMPOSITION` | character | `ADDITIVE`, `DIAGONAL`, `MULTIPLICATIVE`, `SYMMETRIC_MULTIPLICATIVE`, `SCHUR_LOWER`, `SCHUR_UPPER`, `SCHUR_FULL`, `SCHUR_PDE_CONTROL` | global |
 | `BLOCK_SOLVE`, `NEST_BLOCK_SOLVE` | character | `DIAG`, `BJAC`, `NONE` | global or `idx` field |
-| `SCHUR_SOLVE`, `NEST_SCHUR_SOLVE` | character | `A22`, `MATRIX_FREE` | global |
+| `SCHUR_SOLVE`, `NEST_SCHUR_SOLVE` | character | `A22`, `A33`, `MATRIX_FREE`, `SELFP`, `SELF` | global |
 | `SCHUR_MAXIT`, `NEST_SCHUR_MAXIT` | integer | nonnegative iteration count | global |
 | `SCHUR_TOL`, `NEST_SCHUR_TOL` | real | nonnegative tolerance | global |
 | `INNER_SOLVE`, `KRYLOV_SOLVE`, `FIELD_SOLVE` | character | `NONE`, `CG`, `BICGSTAB` | global or `idx` field |
@@ -843,6 +896,12 @@ call prec%set('SCHUR_SOLVE', 'MATRIX_FREE', info)
 call prec%set('SCHUR_MAXIT', 6, info)
 call prec%set('SCHUR_TOL', 1.0d-8, info)
 
+! Alternative two-field Stokes-style Schur approximation:
+! call prec%set('SCHUR_SOLVE', 'SELFP', info)
+! call prec%set('BLOCK_SOLVE', 'BJAC', info)
+! call prec%set('SUB_SOLVE', 'ILU', info)
+! call prec%set('SUB_FILLIN', 0, info)
+
 call prec%set('BLOCK_SOLVE', 'BJAC', info, idx=1)
 call prec%set('BLOCK_SOLVE', 'DIAG', info, idx=2)
 
@@ -857,6 +916,26 @@ call psb_krylov('BICGSTAB', nested_matrix%a_glob, prec, b, x, eps, &
 
 This keeps the outer solve in the existing PSBLAS API while allowing field-specific behavior inside the nested preconditioner.
 
+For a three-field optimal-control KKT system, the PDE-control Schur composition is configured as:
+
+```fortran
+call prec%init('NEST', info)
+call prec%set('COMPOSITION', 'SCHUR_PDE_CONTROL', info)
+call prec%set('SCHUR_SOLVE', 'MATRIX_FREE', info)
+call prec%set('SCHUR_MAXIT', 200, info)
+call prec%set('SCHUR_TOL', 0.0d0, info)
+
+call prec%set('BLOCK_SOLVE', 'DIAG', info, idx=1)
+call prec%set('BLOCK_SOLVE', 'DIAG', info, idx=2)
+call prec%set('BLOCK_SOLVE', 'DIAG', info, idx=3)
+
+call prec%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
+call psb_krylov('BICGSTAB', nested_matrix%a_glob, prec, b, x, eps, &
+     & nested_matrix%desc_glob, info)
+```
+
+The field-3 `DIAG` setting is harmless for zero `A_33` blocks because PSBLAS diagonal preconditioning treats zero diagonal entries as identity. The useful inverse approximation comes from the matrix-free CG solve on the PDE-control Schur action.
+
 ## Factory And API Integration
 
 `psb_dprecinit` recognizes the nested preconditioner name:
@@ -915,12 +994,23 @@ These tests exercise:
 - nested preconditioner initialization
 - additive/diagonal composition
 - multiplicative and symmetric multiplicative composition
-- Schur lower, upper, and full composition
-- matrix-free Schur solve mode
+- Schur lower, upper, full, and PDE-control composition
+- matrix-free Schur solve mode, including CG for `SCHUR_PDE_CONTROL`
 - per-field inner Krylov configuration and apply through the Schur preconditioner
 - serial and MPI execution paths in the nested examples
 
-Validation of the nested preconditioner examples showed convergence for the additive, multiplicative, and Schur-style configurations on the nested Laplacian test problems.
+Validation of the nested preconditioner examples showed convergence for the additive, multiplicative, and Schur-style configurations on the nested Laplacian test problems. The AMG4PSBLAS nested optimal-control MatrixMarket sample also converges with `SCHUR_PDE_CONTROL`, `SCHUR_SOLVE=MATRIX_FREE`, and `SCHUR_MAXIT=200` on serial and two-rank MPI runs.
+
+The PETSc binary Stokes nested test was also used to compare the original matrix-free Schur solve against the new `SELFP` approximation on two MPI ranks. With `SCHUR_FULL`, `BICGSTAB`, `BLOCK_SOLVE=BJAC`, and `SUB_SOLVE=ILU`, the observed runs were:
+
+| Schur solve | Sub fill | Iterations | Relative residual | Solve time |
+| --- | ---: | ---: | ---: | ---: |
+| `MATRIX_FREE` | 1 | 532 | `8.9510E-05` | about `4.95s` |
+| `SELFP` | 0 | 412 | `9.5632E-05` | about `0.50s` |
+| `SELFP` | 1 | 393 | `9.5120E-05` | about `0.51s` |
+| `SELFP` | 2 | 570 | `9.5620E-05` | about `0.80s` |
+
+For this case, `SELFP` reduced the iteration count compared with the matrix-free Schur solve and made each solve substantially cheaper. `SUB_FILLIN=1` gave the lowest iteration count, while `SUB_FILLIN=0` was the best practical default in the measured runs because it was close in iterations and cheapest to build/apply. Larger fill was not automatically better for the lumped Schur approximation.
 
 The Schur nested preconditioner test includes an additional `SCHUR_FULL` / `A22+INNER` case. That case configures independent per-field inner Krylov contexts, builds the nested preconditioner, applies it directly to the right-hand side, and verifies a successful finite preconditioner application. The existing Schur cases continue to check outer BiCGSTAB convergence.
 
@@ -933,8 +1023,9 @@ The current implementation has the following limitations:
 - double precision only
 - nested preconditioner name currently exposed as `NEST`
 - transposed preconditioner application is rejected
-- Schur-style compositions are implemented for exactly two fields
-- matrix-free Schur solve currently uses a small Richardson iteration
+- two-field Schur compositions are implemented by `SCHUR_LOWER`, `SCHUR_UPPER`, and `SCHUR_FULL`; the specialized three-field Schur composition is `SCHUR_PDE_CONTROL` only
+- two-field matrix-free Schur solves use Richardson; PDE-control matrix-free Schur solves use CG on the positive field-3 Schur action
+- `SELFP` is currently a two-field Schur approximation for `SCHUR_LOWER`, `SCHUR_UPPER`, and `SCHUR_FULL`; it assembles a local lumped Schur preconditioning matrix, not an exact distributed Schur complement
 - inner Krylov methods currently include `CG` and `BICGSTAB`
 
 ## Future Work
@@ -942,10 +1033,11 @@ The current implementation has the following limitations:
 Recommended next steps are:
 
 - add single, complex, and complex double precision variants
+- expose a dedicated Schur preconditioner for the PDE-control CG solve, instead of relying on the field-3 block preconditioner fallback
 
 ## Summary
 
-The nested preconditioner adds a field-split preconditioning layer to PSBLAS while preserving the existing solver and preconditioner API. It reuses standard PSBLAS block preconditioners on diagonal field blocks and combines them through additive, multiplicative, symmetric multiplicative, and Schur-style compositions.
+The nested preconditioner adds a field-split preconditioning layer to PSBLAS while preserving the existing solver and preconditioner API. It reuses standard PSBLAS block preconditioners on diagonal field blocks and combines them through additive, multiplicative, symmetric multiplicative, two-field Schur, and three-field PDE-control Schur compositions.
 
 The implementation integrates with:
 

test/nested/README.md

@@ -215,6 +215,33 @@ A field-split interface (`psb_d_nest_get_block`, `psb_d_nest_get_field_desc`, `p
 
 All tests run both serially and in parallel, and the result is invariant with respect to the number of MPI processes.
 
+
+### PETSc Stokes block file test
+
+`psb_d_nest_petsc_stokes_test` reads PETSc-style Stokes MatrixMarket blocks and
+assembles the saddle-point operator as
+
+```text
+[ A   Bt ]
+[ B    0 ]
+```
+
+using `psb_d_nest_matrix`, so the rectangular `Bt` and `B` blocks are kept as
+nested sub-blocks. Configure it with environment variables:
+
+```sh
+export PETSC_STOKES_DIR=/path/to/petsc/stokes/files
+export PETSC_STOKES_METHOD=BICGSTAB
+export PETSC_STOKES_PREC=NEST
+export PETSC_STOKES_COMPOSITION=SCHUR_FULL
+export PETSC_STOKES_SCHUR_SOLVE=MATRIX_FREE
+./runs/psb_d_nest_petsc_stokes_test
+```
+
+The deal.II Stokes file test and optimal-control file test are not part of
+`test/nested`; the retained file-based nested sample here is the PETSc Stokes
+variant.
+
 ### Build and run
 
 The PSBLAS library must be built/installed first (from the repository root):

test/nested/psb_d_nest_petsc_stokes_test.F90

@@ -0,0 +1,446 @@
+!
+! File: psb_d_nest_petsc_stokes_test.F90
+!
+! Read the PETSc binary Stokes files used by PETSc ex87/ex87_stokes,
+! assemble the 2x2 operator as a PSBLAS nested matrix, and solve it through
+! the standard PSBLAS Krylov/preconditioner interface.
+!
+! Expected files in PETSC_STOKES_DIR, or in
+!   ${DATAFILESPATH}/matrices/hpddm/GENEO
+! when PETSC_STOKES_DIR is not set:
+!   B00.dat, B10.dat, B11.dat, rhs_B.dat
+!
+program psb_d_nest_petsc_stokes_test
+  use, intrinsic :: iso_fortran_env, only : int32, real64
+  use psb_base_mod
+  use psb_prec_mod
+  use psb_linsolve_mod
+  use psb_d_nest_mod
+  implicit none
+
+  integer(int32), parameter :: petsc_mat_classid = 1211216_int32
+  integer(int32), parameter :: petsc_vec_classid = 1211214_int32
+
+  type petsc_matrix
+    integer(psb_lpk_) :: nrow = 0_psb_lpk_
+    integer(psb_lpk_) :: ncol = 0_psb_lpk_
+    integer(psb_lpk_) :: nnz  = 0_psb_lpk_
+    integer(psb_lpk_), allocatable :: row(:), col(:)
+    real(psb_dpk_),    allocatable :: val(:)
+  end type petsc_matrix
+
+  type petsc_vector
+    integer(psb_lpk_) :: n = 0_psb_lpk_
+    real(psb_dpk_), allocatable :: val(:)
+  end type petsc_vector
+
+  type(psb_ctxt_type)     :: context
+  type(psb_d_nest_matrix) :: nested_matrix
+  type(psb_dprec_type)    :: preconditioner
+  type(psb_d_vect_type)   :: rhs, x_solution, residual
+
+  type(petsc_matrix) :: b00, b10, b11
+  type(petsc_vector) :: rhs_full
+
+  integer(psb_ipk_) :: my_rank, num_procs, info
+  integer(psb_ipk_) :: max_iter, trace_level, stop_criterion, n_iter, n_insert, schur_maxit, sub_fillin
+  integer(psb_lpk_) :: n_u, n_p, n_total
+  integer(psb_lpk_), allocatable :: rows(:), cols(:)
+  real(psb_dpk_),    allocatable :: vals(:)
+  real(psb_dpk_) :: eps, check_tol, final_residual, residual_norm, rhs_norm
+  real(psb_dpk_) :: t0, t_read, t_assemble, t_prec, t_solve
+  character(len=256) :: input_dir, method, ptype, composition, schur_solve, block_solve, sub_solve
+  logical :: empty_a11
+
+  call psb_init(context)
+  call psb_info(context, my_rank, num_procs)
+
+  call get_stokes_dir(input_dir)
+  call get_string_env('PETSC_STOKES_METHOD', method, 'BICGSTAB')
+  call get_string_env('PETSC_STOKES_PREC', ptype, 'NEST')
+  call get_string_env('PETSC_STOKES_COMPOSITION', composition, 'SCHUR_FULL')
+  call get_string_env('PETSC_STOKES_SCHUR_SOLVE', schur_solve, 'MATRIX_FREE')
+  call get_string_env('PETSC_STOKES_BLOCK_SOLVE', block_solve, 'BJAC')
+  call get_string_env('PETSC_STOKES_SUB_SOLVE', sub_solve, 'ILU')
+  call get_int_env('PETSC_STOKES_ITMAX', max_iter, 2000)
+  call get_int_env('PETSC_STOKES_TRACE', trace_level, 0)
+  call get_int_env('PETSC_STOKES_SCHUR_MAXIT', schur_maxit,20)
+  call get_int_env('PETSC_STOKES_SUB_FILLIN', sub_fillin, 1)
+  call get_real_env('PETSC_STOKES_EPS', eps, 1.0e-4_psb_dpk_)
+  call get_real_env('PETSC_STOKES_CHECK_TOL', check_tol, 1.0e-4_psb_dpk_)
+  stop_criterion = 2
+
+  if (my_rank == 0) then
+    write(*,'(a)') 'PETSc Stokes nested binary test'
+    write(*,'(a,a)') '  input dir  : ', trim(input_dir)
+    write(*,'(a,a)') '  method     : ', trim(method)
+    write(*,'(a,a)') '  prec       : ', trim(ptype)
+    write(*,'(a,a)') '  composition: ', trim(composition)
+    write(*,'(a,a)') '  schur solve: ', trim(schur_solve)
+    write(*,'(a,a)') '  block solve: ', trim(block_solve)
+    write(*,'(a,a)') '  sub solve  : ', trim(sub_solve)
+    write(*,'(a,i0)') '  schur maxit: ', schur_maxit
+    write(*,'(a,i0)') '  sub fillin : ', sub_fillin
+    write(*,'(a,es12.4)') '  tolerance  : ', eps
+    write(*,'(a,es12.4)') '  check tol  : ', check_tol
+  end if
+
+  t0 = psb_wtime()
+  call read_petsc_matrix(path_join(input_dir, 'B00.dat'), b00)
+  call read_petsc_matrix(path_join(input_dir, 'B10.dat'), b10)
+  call read_petsc_matrix(path_join(input_dir, 'B11.dat'), b11)
+  call read_petsc_vector(path_join(input_dir, 'rhs_B.dat'), rhs_full)
+  t_read = psb_wtime() - t0
+
+  n_u = b00%nrow
+  n_p = b10%nrow
+  n_total = n_u + n_p
+  empty_a11 = (b11%nnz == 0_psb_lpk_)
+
+  if ((b00%ncol /= n_u) .or. (b10%ncol /= n_u) .or. &
+      (b11%nrow /= n_p) .or. (b11%ncol /= n_p) .or. &
+      (rhs_full%n /= n_total)) then
+    if (my_rank == 0) write(*,*) 'FAIL: inconsistent PETSc Stokes dimensions'
+    call psb_abort(context)
+  end if
+
+  t0 = psb_wtime()
+  call nested_matrix%init(context, [n_u, n_p], info)
+  call check_info(info, 'nested_matrix%init')
+
+  call select_owned_entries(b00, nested_matrix%get_owned_rows(1), rows, cols, vals, n_insert)
+  call nested_matrix%ins(1, 1, n_insert, rows, cols, vals, info)
+  call check_info(info, 'insert B00')
+  call clear_triplets(rows, cols, vals)
+
+  call select_owned_entries_transpose(b10, nested_matrix%get_owned_rows(1), rows, cols, vals, n_insert)
+  call nested_matrix%ins(1, 2, n_insert, rows, cols, vals, info)
+  call check_info(info, 'insert B10 transpose')
+  call clear_triplets(rows, cols, vals)
+
+  call select_owned_entries(b10, nested_matrix%get_owned_rows(2), rows, cols, vals, n_insert)
+  call nested_matrix%ins(2, 1, n_insert, rows, cols, vals, info)
+  call check_info(info, 'insert B10')
+  call clear_triplets(rows, cols, vals)
+
+  if (.not. empty_a11) then
+    call select_owned_entries(b11, nested_matrix%get_owned_rows(2), rows, cols, vals, n_insert)
+    call nested_matrix%ins(2, 2, n_insert, rows, cols, vals, info)
+    call check_info(info, 'insert B11')
+    call clear_triplets(rows, cols, vals)
+  end if
+
+  call nested_matrix%asb(info)
+  call check_info(info, 'nested_matrix%asb')
+  t_assemble = psb_wtime() - t0
+
+  call psb_geall(rhs, nested_matrix%desc_glob, info)
+  call check_info(info, 'psb_geall(rhs)')
+  call insert_owned_rhs(rhs_full, n_u, nested_matrix%get_owned_rows(1), nested_matrix%get_owned_rows(2), &
+       & rhs, nested_matrix%desc_glob, info)
+  call check_info(info, 'insert rhs')
+  call psb_geasb(rhs, nested_matrix%desc_glob, info, dupl=psb_dupl_add_)
+  call check_info(info, 'psb_geasb(rhs)')
+
+  call psb_geall(x_solution, nested_matrix%desc_glob, info)
+  call check_info(info, 'psb_geall(x)')
+  call x_solution%zero()
+  call psb_geasb(x_solution, nested_matrix%desc_glob, info)
+  call check_info(info, 'psb_geasb(x)')
+
+  t0 = psb_wtime()
+  call preconditioner%init(context, trim(ptype), info)
+  call check_info(info, 'prec%init')
+  if (psb_toupper(trim(ptype)) == 'NEST') then
+    call preconditioner%set('COMPOSITION', trim(composition), info)
+    call preconditioner%set('SCHUR_SOLVE', trim(schur_solve), info)
+    call preconditioner%set('SCHUR_MAXIT', schur_maxit, info)
+    call preconditioner%set('BLOCK_SOLVE', trim(block_solve), info)
+    call preconditioner%set('SUB_SOLVE', trim(sub_solve), info)
+    call preconditioner%set('SUB_FILLIN', sub_fillin, info)
+  end if
+  call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
+  call check_info(info, 'prec%build')
+  t_prec = psb_wtime() - t0
+
+  t0 = psb_wtime()
+  call psb_krylov(trim(method), nested_matrix%a_glob, preconditioner, rhs, x_solution, eps, &
+       & 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')
+  t_solve = psb_wtime() - t0
+
+  call psb_geall(residual, nested_matrix%desc_glob, info)
+  call psb_geasb(residual, nested_matrix%desc_glob, info)
+  call psb_geaxpby(done, rhs, dzero, residual, nested_matrix%desc_glob, info)
+  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)
+  rhs_norm      = psb_genrm2(rhs,      nested_matrix%desc_glob, info)
+
+  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
+    write(*,'(a,l1)') '  empty B11  : ', empty_a11
+    write(*,'(a,i8,a,es12.4)') '  iterations=', n_iter, '  solver err=', final_residual
+    write(*,'(a,es12.4)') '  ||b-Ax||_2 / ||b||_2=', residual_norm / max(rhs_norm, tiny(done))
+    if (residual_norm / max(rhs_norm, tiny(done)) > check_tol) then
+      write(*,'(a)') '[FAIL] psb_d_nest_petsc_stokes_test: residual above PETSC_STOKES_CHECK_TOL'
+      call psb_abort(context)
+    end if
+    write(*,'(a,es12.4,a,es12.4,a,es12.4,a,es12.4)') &
+         & '  times read=', t_read, ' assemble=', t_assemble, ' prec=', t_prec, ' solve=', t_solve
+    write(*,'(a)') '[PASS] psb_d_nest_petsc_stokes_test'
+  end if
+
+9999 continue
+  call preconditioner%free(info)
+  call psb_gefree(residual, nested_matrix%desc_glob, info)
+  call psb_gefree(x_solution, nested_matrix%desc_glob, info)
+  call psb_gefree(rhs, nested_matrix%desc_glob, info)
+  call nested_matrix%free(info)
+  call psb_exit(context)
+
+contains
+
+  subroutine check_info(info_value, label)
+    integer(psb_ipk_), intent(in) :: info_value
+    character(len=*), intent(in) :: label
+    if (info_value /= psb_success_) then
+      if (my_rank == 0) write(*,*) 'FAIL: ', trim(label), ' info=', info_value
+      call psb_abort(context)
+    end if
+  end subroutine check_info
+
+  function path_join(dir, name) result(path)
+    character(len=*), intent(in) :: dir, name
+    character(len=512) :: path
+    integer :: n
+    n = len_trim(dir)
+    if (n > 0 .and. dir(n:n) == '/') then
+      path = trim(dir) // trim(name)
+    else
+      path = trim(dir) // '/' // trim(name)
+    end if
+  end function path_join
+
+  subroutine get_stokes_dir(value)
+    character(len=*), intent(out) :: value
+    character(len=256) :: datafiles
+    integer :: status
+
+    call get_environment_variable('PETSC_STOKES_DIR', value, status=status)
+    if (status == 0 .and. len_trim(value) > 0) return
+
+    call get_environment_variable('DATAFILESPATH', datafiles, status=status)
+    if (status == 0 .and. len_trim(datafiles) > 0) then
+      value = trim(datafiles) // '/matrices/hpddm/GENEO'
+    else
+      value = '.'
+    end if
+  end subroutine get_stokes_dir
+
+  subroutine get_string_env(name, value, default_value)
+    character(len=*), intent(in) :: name, default_value
+    character(len=*), intent(out) :: value
+    integer :: status
+    call get_environment_variable(name, value, status=status)
+    if (status /= 0 .or. len_trim(value) == 0) value = default_value
+  end subroutine get_string_env
+
+  subroutine get_int_env(name, value, default_value)
+    character(len=*), intent(in) :: name
+    integer(psb_ipk_), intent(out) :: value
+    integer(psb_ipk_), intent(in) :: default_value
+    character(len=64) :: text
+    integer :: status
+    call get_environment_variable(name, text, status=status)
+    if (status == 0 .and. len_trim(text) > 0) then
+      read(text,*) value
+    else
+      value = default_value
+    end if
+  end subroutine get_int_env
+
+  subroutine get_real_env(name, value, default_value)
+    character(len=*), intent(in) :: name
+    real(psb_dpk_), intent(out) :: value
+    real(psb_dpk_), intent(in) :: default_value
+    character(len=64) :: text
+    integer :: status
+    call get_environment_variable(name, text, status=status)
+    if (status == 0 .and. len_trim(text) > 0) then
+      read(text,*) value
+    else
+      value = default_value
+    end if
+  end subroutine get_real_env
+
+  subroutine read_petsc_matrix(filename, mat)
+    character(len=*), intent(in) :: filename
+    type(petsc_matrix), intent(out) :: mat
+    integer :: unit, ios
+    integer(int32) :: classid, nrow32, ncol32, nnz32
+    integer(int32), allocatable :: row_nnz(:), cols0(:)
+    real(real64), allocatable :: vals64(:)
+    integer(psb_lpk_) :: i, j, k
+
+    open(newunit=unit, file=trim(filename), status='old', action='read', access='stream', &
+         & form='unformatted', convert='big_endian', iostat=ios)
+    if (ios /= 0) then
+      write(*,*) 'Could not open PETSc matrix file: ', trim(filename)
+      stop 2
+    end if
+    read(unit) classid
+    if (classid /= petsc_mat_classid) stop 'Bad PETSc matrix classid'
+    read(unit) nrow32
+    read(unit) ncol32
+    read(unit) nnz32
+    mat%nrow = int(nrow32, psb_lpk_)
+    mat%ncol = int(ncol32, psb_lpk_)
+    mat%nnz  = int(nnz32,  psb_lpk_)
+
+    allocate(row_nnz(mat%nrow))
+    if (mat%nrow > 0_psb_lpk_) read(unit) row_nnz
+    allocate(cols0(mat%nnz), vals64(mat%nnz))
+    if (mat%nnz > 0_psb_lpk_) then
+      read(unit) cols0
+      read(unit) vals64
+    end if
+    close(unit)
+
+    allocate(mat%row(mat%nnz), mat%col(mat%nnz), mat%val(mat%nnz))
+    k = 0_psb_lpk_
+    do i = 1_psb_lpk_, mat%nrow
+      do j = 1_psb_lpk_, int(row_nnz(i), psb_lpk_)
+        k = k + 1_psb_lpk_
+        mat%row(k) = i
+        mat%col(k) = int(cols0(k), psb_lpk_) + 1_psb_lpk_
+        mat%val(k) = real(vals64(k), psb_dpk_)
+      end do
+    end do
+    if (k /= mat%nnz) stop 'Bad PETSc matrix row lengths'
+  end subroutine read_petsc_matrix
+
+  subroutine read_petsc_vector(filename, vec)
+    character(len=*), intent(in) :: filename
+    type(petsc_vector), intent(out) :: vec
+    integer :: unit, ios
+    integer(int32) :: classid, n32
+    real(real64), allocatable :: vals64(:)
+
+    open(newunit=unit, file=trim(filename), status='old', action='read', access='stream', &
+         & form='unformatted', convert='big_endian', iostat=ios)
+    if (ios /= 0) then
+      write(*,*) 'Could not open PETSc vector file: ', trim(filename)
+      stop 2
+    end if
+    read(unit) classid
+    if (classid /= petsc_vec_classid) stop 'Bad PETSc vector classid'
+    read(unit) n32
+    vec%n = int(n32, psb_lpk_)
+    allocate(vals64(vec%n), vec%val(vec%n))
+    if (vec%n > 0_psb_lpk_) read(unit) vals64
+    close(unit)
+    vec%val = real(vals64, psb_dpk_)
+  end subroutine read_petsc_vector
+
+  subroutine select_owned_entries(mat, owned_rows, out_rows, out_cols, out_vals, count)
+    type(petsc_matrix), intent(in) :: mat
+    integer(psb_lpk_), intent(in) :: owned_rows(:)
+    integer(psb_lpk_), allocatable, intent(out) :: out_rows(:), out_cols(:)
+    real(psb_dpk_), allocatable, intent(out) :: out_vals(:)
+    integer(psb_ipk_), intent(out) :: count
+    integer(psb_lpk_) :: i, first_owned, last_owned
+    integer(psb_ipk_) :: k
+
+    count = 0
+    if (size(owned_rows) == 0) then
+      allocate(out_rows(0), out_cols(0), out_vals(0))
+      return
+    end if
+
+    first_owned = minval(owned_rows)
+    last_owned  = maxval(owned_rows)
+    do i = 1_psb_lpk_, mat%nnz
+      if (mat%row(i) >= first_owned .and. mat%row(i) <= last_owned) count = count + 1
+    end do
+
+    allocate(out_rows(count), out_cols(count), out_vals(count))
+    k = 0
+    do i = 1_psb_lpk_, mat%nnz
+      if (mat%row(i) >= first_owned .and. mat%row(i) <= last_owned) then
+        k = k + 1
+        out_rows(k) = mat%row(i)
+        out_cols(k) = mat%col(i)
+        out_vals(k) = mat%val(i)
+      end if
+    end do
+  end subroutine select_owned_entries
+
+  subroutine select_owned_entries_transpose(mat, owned_rows, out_rows, out_cols, out_vals, count)
+    type(petsc_matrix), intent(in) :: mat
+    integer(psb_lpk_), intent(in) :: owned_rows(:)
+    integer(psb_lpk_), allocatable, intent(out) :: out_rows(:), out_cols(:)
+    real(psb_dpk_), allocatable, intent(out) :: out_vals(:)
+    integer(psb_ipk_), intent(out) :: count
+    integer(psb_lpk_) :: i, first_owned, last_owned
+    integer(psb_ipk_) :: k
+
+    count = 0
+    if (size(owned_rows) == 0) then
+      allocate(out_rows(0), out_cols(0), out_vals(0))
+      return
+    end if
+
+    first_owned = minval(owned_rows)
+    last_owned  = maxval(owned_rows)
+    do i = 1_psb_lpk_, mat%nnz
+      if (mat%col(i) >= first_owned .and. mat%col(i) <= last_owned) count = count + 1
+    end do
+
+    allocate(out_rows(count), out_cols(count), out_vals(count))
+    k = 0
+    do i = 1_psb_lpk_, mat%nnz
+      if (mat%col(i) >= first_owned .and. mat%col(i) <= last_owned) then
+        k = k + 1
+        out_rows(k) = mat%col(i)
+        out_cols(k) = mat%row(i)
+        out_vals(k) = mat%val(i)
+      end if
+    end do
+  end subroutine select_owned_entries_transpose
+
+  subroutine insert_owned_rhs(vec, velocity_size, owned_u, owned_p, rhs_vec, desc, info)
+    type(petsc_vector), intent(in) :: vec
+    integer(psb_lpk_), intent(in) :: velocity_size
+    integer(psb_lpk_), intent(in) :: owned_u(:), owned_p(:)
+    type(psb_d_vect_type), intent(inout) :: rhs_vec
+    type(psb_desc_type), intent(in) :: desc
+    integer(psb_ipk_), intent(out) :: info
+    integer(psb_lpk_), allocatable :: idx(:)
+    real(psb_dpk_), allocatable :: values(:)
+    integer :: i, n
+
+    n = size(owned_u) + size(owned_p)
+    allocate(idx(n), values(n))
+    do i = 1, size(owned_u)
+      idx(i) = owned_u(i)
+      values(i) = vec%val(owned_u(i))
+    end do
+    do i = 1, size(owned_p)
+      idx(size(owned_u) + i) = velocity_size + owned_p(i)
+      values(size(owned_u) + i) = vec%val(velocity_size + owned_p(i))
+    end do
+    call psb_geins(n, idx, values, rhs_vec, desc, info)
+    deallocate(idx, values)
+  end subroutine insert_owned_rhs
+
+  subroutine clear_triplets(in_rows, in_cols, in_vals)
+    integer(psb_lpk_), allocatable, intent(inout) :: in_rows(:), in_cols(:)
+    real(psb_dpk_), allocatable, intent(inout) :: in_vals(:)
+    if (allocated(in_rows)) deallocate(in_rows)
+    if (allocated(in_cols)) deallocate(in_cols)
+    if (allocated(in_vals)) deallocate(in_vals)
+  end subroutine clear_triplets
+
+end program psb_d_nest_petsc_stokes_test

test/nested/psb_d_nest_schur_prec_cg_test.F90

@@ -302,41 +302,12 @@ program psb_d_nest_schur_prec_cg_test
   call psb_gefree(x_solution, nested_matrix%desc_glob, info)
   call preconditioner%free(info)
 
-  ! Exercise transposed Schur application for the explicit A22 approximation.
-  call preconditioner%init(context, 'NEST', info)
-  call preconditioner%set('COMPOSITION', 'SCHUR_FULL', info)
-  call preconditioner%set('SCHUR_SOLVE', 'A22', info)
-  call preconditioner%set('BLOCK_SOLVE', 'DIAG', info)
-  call preconditioner%build(nested_matrix%a_glob, nested_matrix%desc_glob, info)
-  if (info /= psb_success_) then
-    if (my_rank == 0) write(*,*) 'FAIL: nested Schur transpose prec%build info=', info
-    all_passed = .false.
-    goto 9999
-  end if
-
-  call psb_geall(x_solution, nested_matrix%desc_glob, info)
-  call psb_geasb(x_solution, nested_matrix%desc_glob, info)
-  call preconditioner%apply(rhs, x_solution, nested_matrix%desc_glob, info, trans='T')
-  if (info /= psb_success_) then
-    if (my_rank == 0) write(*,*) 'FAIL: NEST Schur transpose apply info=', info
-    all_passed = .false.
-  end if
-
-  solution_error = psb_genrm2(x_solution, nested_matrix%desc_glob, info)
+  ! Transposed nested preconditioner application is intentionally unsupported.
+  ! The implementation returns psb_err_transpose_not_n_unsupported_ for trans /= 'N'.
   if (my_rank == 0) then
-    write(*,'(a,a32,a,a16,a,a16,a,es12.4)') &
-         & ' prec=NEST/', 'SCHUR_FULL', &
-         & ' schur=', 'A22', &
-         & ' method=', 'APPLY^T', &
-         & '  ||P^T rhs||=', solution_error
+    write(*,'(a)') ' prec=NEST/SCHUR_FULL schur=A22 method=APPLY^T skipped: transpose unsupported'
   end if
 
-  if ((info /= psb_success_) .or. (solution_error /= solution_error) .or. &
-       & (solution_error <= dzero)) all_passed = .false.
-
-  call psb_gefree(x_solution, nested_matrix%desc_glob, info)
-  call preconditioner%free(info)
-
   if (my_rank == 0) then
     if (all_passed) then
       write(*,*) '[PASS] Krylov solvers converge with Schur-style NEST preconditioners'