Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse sparse_jac sparse_jac_xam.m

cppad_swig-> testing swig_xam.i library whats_new_2018

library-> module a_double vector a_fun sparse error xam.m4

sparse-> sparse_rc sparse_rcv jac_sparsity hes_sparsity sparse_jac sparse_hes

sparse_jac-> sparse_jac_xam.cpp sparse_jac_xam.m sparse_jac_xam.pm sparse_jac_xam.py

sparse_jac_xam.m

Headings

function ok = sparse_jac_xam()
     %
     % load the Cppad Swig library
     m_cppad
     %
     % initialize return variable
     ok = true;
     % -----------------------------------------------------------------------
     % number of dependent and independent variables
     n = 3;
     % one
     aone = m_cppad.a_double(1.0);
     %
     % create the independent variables ax
     x = m_cppad.vec_double(n);
     for i = [ 0 :(n -1) ]
          x(i) = i + 2.0;
     end
     ax = m_cppad.independent(x);
     %
     % create dependent variables ay with ay[i] = (j+1) * ax[j]
     % where i = mod(j + 1, n)
     ay = m_cppad.vec_a_double(n);
     for j = [ 0 :(n -1) ]
          i = j+1;
          if( i >= n  )
               i = i - n;
          end
          aj = m_cppad.a_double(j);
          ay_i = (aj + aone) * ax(j);
          ay(i) = ay_i;
     end
     %
     % define af corresponding to f(x)
     af = m_cppad.a_fun(ax, ay);
     %
     % sparsity pattern for identity matrix
     pat_eye = m_cppad.sparse_rc();
     pat_eye.resize(n, n, n);
     for k = [ 0 :(n-1) ]
          pat_eye.put(k, k, k);
     end
     %
     % sparsity pattern for the Jacobian
     pat_jac = m_cppad.sparse_rc();
     af.for_jac_sparsity(pat_eye, pat_jac);
     %
     % loop over forward and reverse mode
     for mode = [ 0 :(2-1) ]
          % compute all possibly non-zero entries in Jacobian
          subset = m_cppad.sparse_rcv(pat_jac);
          % work space used to save time for multiple calls
          work = m_cppad.sparse_jac_work();
          if( mode == 0  )
               af.sparse_jac_for(subset, x, pat_jac, work);
          end
          if( mode == 1  )
               af.sparse_jac_rev(subset, x, pat_jac, work);
          end
          %
          % check result
          ok = ok && n == subset.nnz();
          col_major = subset.col_major();
          row = subset.row();
          col = subset.col();
          val = subset.val();
          for k = [ 0 :(n-1) ]
               ell = col_major(k);
               r = row(ell);
               c = col(ell);
               v = val(ell);
               i = c+1;
               if( i >=  n  )
                    i = i - n;
               end
               ok = ok && c == k;
               ok = ok && r == i;
               ok = ok && v == c + 1.0;
          end
     end
     %
     return;
end