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.