Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse sparse_jac sparse_jac_xam.py

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.py

Headings

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