-------------------------------------------------------------------------------
lines 8-179 of file: include/cppad/example/atomic_four/lin_ode/jac_sparsity.hpp
-------------------------------------------------------------------------------

{xrst_begin atomic_four_lin_ode_jac_sparsity.hpp}
{xrst_spell
   nnz
}

Atomic Linear ODE Jacobian Sparsity Pattern: Example Implementation
###################################################################

Purpose
*******
The ``jac_sparsity`` routine overrides the virtual functions
used by the atomic_four base class for Jacobian sparsity calculations; see
:ref:`jac_sparsity<atomic_four_jac_sparsity-name>` .

Notation
********
We use the notation:
:ref:`atomic_four_lin_ode@call_id`
:ref:`atomic_four_lin_ode@r`
:ref:`atomic_four_lin_ode@pattern`
:ref:`atomic_four_lin_ode@transpose`
:ref:`atomic_four_lin_ode@pattern@nnz` ,
:ref:`atomic_four_lin_ode@pattern@row` ,
:ref:`atomic_four_lin_ode@pattern@col` ,
:ref:`atomic_four_lin_ode@x` ,
:ref:`atomic_four_lin_ode@x@n` ,
:ref:`atomic_four_lin_ode@x@A(x)` ,
:ref:`atomic_four_lin_ode@x@b(x)` ,
:ref:`atomic_four_lin_ode@y(x)` ,
:ref:`atomic_four_lin_ode@y(x)@m` ,
:ref:`atomic_four_lin_ode@vk(x)` ,
and the following additional notation:

S[ g(x) ]
=========
We use :math:`S [ g(x) ]` to denote the sparsity pattern
for a function :math:`g : \B{R}^n \rightarrow \B{R}^m` as a vector of sets.
To be specific, for :math:`i = 0, \ldots , m-1`,
:math:`S_i [ g(x) ]` is the set of indices between
zero and :math:`n - 1` such that
:math:`\partial g_i (x) / \partial x_j` is possibly non-zero.

N [ g(x) ]
==========
We use :math:`N[ g(x) ]` to denote the set of :math:`i`
such that :math:`g_i (x)` is not identically zero.

J_i [ A(x) ]
============
We use :math:`J_i [ A(x) ]` to denote the set of :math:`j`
between zero and :math:`m-1` such that
:math:`A_{i,j}` is not known to be identically zero.

P_i [ g(x) ]
============
We use :math:`P_i [ g(x) ]` to denote the set if sparsity pattern indices

.. math::

   P_i [ g(x) ] = \left\{ p  \W{:}
      0 \leq p < \R{nnz} \W{,}
      i = \R{row} [p] \W{,}
      \R{col}[p] \in N [ g(x) ]
   \right\}

Theory
******
This routine must calculate the following value for
:math:`i = 0, \ldots, m-1`; see :ref:`atomic_four_lin_ode@y(x)@m` :

.. math::

   S_i [ y (x) ] = \bigcup_k S_i [  v^k (x) ]

The set :math:`S_i [ v^0 (x) ]` has just one element
corresponding to :math:`b_i (x)`; i.e,

.. math::

   S_i [ v^0 (x) ] = \{ \R{nnz} + i \}

see :ref:`atomic_four_lin_ode@x@b(x)` .
Furthermore, for :math:`k > 0`,

.. math::

   v^k (x) = \frac{r}{k} A(x) v^{k-1} (x)

.. math::

   S_i [ v^k (x) ] = S_i [ A(x) v^{k-1} (x) ]

.. math::

   S_i [ v^k (x) ] =  P_i [ v^{k-1} (x) ]
   \cup \left\{ S_j [ v^{k-1} (x) ] \W{:}  j \in J_i [ A(x) ] \right\}

Suppose that :math:`\ell` is such that for all :math:`i`
the following two conditions hold

.. math::

   N [ v^\ell (x) ]  \subset \bigcup_{k < \ell} N [ v^k (x) ]

.. math::

   S_i [ v^\ell (x) ] \subset \bigcup_{k < \ell} S_i [ v^k (x) ]

From the first condition above it follows that

.. math::

   P_i [ v^\ell (x) ] \subset \bigcup_{k < \ell} P_i [ v^k (x) ]

Using the second condition we have

.. math::

   S_i [ v^{\ell+1} (x) ] =  P_i [ v^\ell (x) ]
   \cup \left\{ S_j [ v^\ell (x) ] \W{:} j \in J_i [ A(x) ] \right\}

.. math::

   S_i [ v^{\ell+1} (x) ] \subset
   \left\{ S_j [ v^\ell (x) ] \W{:} j \in J_i [ A(x) ] \right\}
      \bigcup_{k \leq \ell } S_i [ v^k (x) ]

.. math::

   S_i [ v^{\ell+1} (x) ]
   \subset
   \bigcup_{k \leq \ell} S_i [ v^k (x) ] \cup
         \left\{ S_j [ v^k (x) ] \W{:} j \in J_i [ A(x) ] \right\}

.. math::

   S_i [ v^{\ell+1} (x) ]
   \subset
   \bigcup_{k < \ell} S_i [ v^k (x) ] \cup
         \left\{ S_j [ v^k (x) ] \W{:} j \in J_i [ A(x) ] \right\}

.. math::

   S_i [ v^{\ell+1} (x) ]
   \subset
   \bigcup_{k \leq \ell} S_i [ v^k (x) ]

.. math::

   S_i [ v^{\ell+1} (x) ]
   \subset
   \bigcup_{k < \ell} S_i [ v^k (x) ]

It follows that

.. math::

   S_i [ y(x) ] = \bigcup_{k < \ell} S_i [ v^k (x) ]

Example
*******
The file :ref:`atomic_four_lin_ode_sparsity.cpp-name`
contains an example and test using this operator.

Source
******
{xrst_literal
   // BEGIN C++
   // END C++
}

{xrst_end atomic_four_lin_ode_jac_sparsity.hpp}