I'm looking for a method to solve a 20000 variable, 20000 residual non-linear PDE with a Galerkin method.
I have Fortran subroutines for:
The residuals: $\vec{r}(\vec{x})$;
Their Jacobian multiplied by a direct algorithmic derivative seed: $\dot{r}(\dot{x})=\frac{\partial \vec{r}(\vec{x})}{\partial \vec{x}}\cdot\dot{x}$; and
Their transposed Jacobian multiplied by a reverse algorithmic derivative seed:
$$\bar{x}(\bar{r})=\left[\frac{\partial \vec{r}(\vec{x})}{\partial \vec{x}}\right]^t\cdot \bar{r}\, .$$
All these subroutines are run in no more than 0.05 s each. Due to the high number of variables, however, I cannot fully compute the Jacobian for a Newton-Raphson iteration because the run time would become impractical.
Is there any out-of-the-shelf method I can use in this case to solve the equation without fully calculating the Jacobian?