<jupyter_start><jupyter_text>Restarting Linear Solutions for Expensive Linear SolvesWhen using iterative linear solvers, it is often desirable to use the converged solution from a previous linear solve as the initial guess for the current one. There is some memory cost associated with this feature, because the solution for each quantity of interest will be saved separately. However, the benefit is reduced computational cost for the subsequent linear solves.```{Note}This feature should not be used when using the [DirectSolver](../../building_blocks/solvers/direct_solver) at the top level of your model. It won’t offer any computational savings in that situation.```To use this feature, provide `cache_linear_solution=True` as an argument to [add_design_var()](../adding_desvars_cons_objs/adding_design_variables), [add_objective()](../adding_desvars_cons_objs/adding_objective) , or [add_constraint()](../adding_desvars_cons_objs/adding_constraint).If you are using one of the OpenMDAO iterative solvers ([ScipyKrylov](../../building_blocks/solvers/scipy_iter_solver), [PETScKrylov](../../building_blocks/solvers/petsc_krylov), [LinearBlockGS](../../building_blocks/solvers/linear_block_gs), or [LinearBlockJac](../../building_blocks/solvers/linear_block_jac), then simply adding that argument is enough to activate this feature.If you have implemented the `solve_linear`() method for an [ImplicitComponent](../working_with_components/implicit_component), then you will need to make sure to use the provided guess solution in your implementation. The cached solution will be put into the solution vector for you to use as an initial guess. Note that you will override those values with the final solution. In `fwd` mode, the guess will be in the `d_outputs` vector. In `rev` mode, the guess will be in the `d_residuals` vector.Below is a toy example problem that illustrates how the restart vectors should work. The restart is passed in via the `x0` argument to gmres.<jupyter_code>from distutils.version import LooseVersion
import numpy as np
import scipy
from scipy.sparse.linalg import gmres

import openmdao.api as om


class QuadraticComp(om.ImplicitComponent):
    """
    A Simple Implicit Component representing a Quadratic Equation.

    R(a, b, c, x) = ax^2 + bx + c

    Solution via Quadratic Formula:
    x = (-b + sqrt(b^2 - 4ac)) / 2a
    """

    def setup(self):
        self.add_input('a', val=1.)
        self.add_input('b', val=1.)
        self.add_input('c', val=1.)
        self.add_output('states', val=[0,0])

    def setup_partials(self):
        self.declare_partials(of='*', wrt='*')

    def apply_nonlinear(self, inputs, outputs, residuals):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['states'][0]
        y = outputs['states'][1]

        residuals['states'][0] = a * x ** 2 + b * x + c
        residuals['states'][1] = a * y + b

    def solve_nonlinear(self, inputs, outputs):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        outputs['states'][0] = (-b + (b ** 2 - 4 * a * c) ** 0.5) / (2 * a)
        outputs['states'][1] = -b/a

    def linearize(self, inputs, outputs, partials):
        a = inputs['a'][0]
        b = inputs['b'][0]
        c = inputs['c'][0]
        x = outputs['states'][0]
        y = outputs['states'][1]

        partials['states', 'a'] = [[x**2],[y]]
        partials['states', 'b'] = [[x],[1]]
        partials['states', 'c'] = [[1.0],[0]]
        partials['states', 'states'] = [[2*a*x+b, 0],[0, a]]

        self.state_jac = np.array([[2*a*x+b, 0],[0, a]])

    def solve_linear(self, d_outputs, d_residuals, mode):

        if mode == 'fwd':
            print("incoming initial guess", d_outputs['states'])
            d_outputs['states'] = gmres(self.state_jac, d_residuals['states'], x0=d_outputs['states'])[0]

        elif mode == 'rev':
            d_residuals['states'] = gmres(self.state_jac, d_outputs['states'], x0=d_residuals['states'])[0]

p = om.Problem()
p.driver = om.ScipyOptimizeDriver()
p.driver.options['optimizer'] = 'SLSQP'

p.model.set_input_defaults('b', val=4.)
p.model.add_subsystem('quad', QuadraticComp(), promotes_inputs=['a', 'b', 'c'], promotes_outputs=['states'])
p.model.add_subsystem('obj', om.ExecComp('y = (x[1]-x[0])**2', x=np.ones(2)))
p.model.connect('states', 'obj.x')

p.model.add_design_var('a', upper=4., cache_linear_solution=True)
p.model.add_constraint('states', upper=10)
p.model.add_objective('obj.y')

p.setup(mode='fwd')
p.run_driver()

print(p['a'], p['b'], p['c'])
print(p['states'])
print(p['obj.y'])
from openmdao.utils.assert_utils import assert_near_equal
import numpy as np

assert_near_equal(p['a'], 3.99999858, tolerance=1e-3)
assert_near_equal(p['b'], 4.0, tolerance=1e-3)
assert_near_equal(p['c'], 1.0, tolerance=1e-3)
assert_near_equal(p['states'], np.array([-0.49970278, -1.00000035]), tolerance=1e-3)
assert_near_equal(p['obj.y'], 0.25000053, tolerance=1e-3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15503100620124025, 'fraction_numerical': 0.020004000800160033, 'mean_word_length': 4.08130081300813, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '389082', 'license': 'Apache-2.0', 'n_tokens_mistral': 1719, 'n_tokens_neox': 1642, 'n_words': 447, 'path': 'openmdao/docs/openmdao_book/features/core_features/working_with_derivatives/linear_restart.ipynb', 'repo_name': 'tong0711/OpenMDAO', 'size': '4999'}