  We consider parametric families of partial differential equations--PDEs where
the parameter $\kappa$ modifies only the (1,1) block of a saddle point matrix
product of a discretization below. The main goal is to develop an algorithm
that removes, as much as possible, the dependence of iterative solvers on the
parameter $\kappa$. The algorithm we propose requires only one matrix
factorization which does not depend on $\kappa$, therefore, allows to reuse it
for solving very fast a large number of discrete PDEs for different $\kappa$
and forcing terms. The design of the proposed algorithm is motivated by
previous works on natural factor of formulation of the stiffness matrices and
their stable numerical solvers. As an application, in two dimensions, we
consider an iterative preconditioned solver based on the null space of
Crouzeix-Raviart discrete gradient represented as the discrete curl of $P_1$
conforming finite element functions. For the numerical examples, we consider
the case of random coefficient pressure equation where the permeability is
modeled by an stochastic process. We note that contrarily from recycling Krylov
subspace techniques, the proposed algorithm does not require fixed forcing
terms.
