|
|
Prev | Next | sparse_hes_xam.py | Headings |
def sparse_hes_xam() :
#
# load the Cppad Swig library
import py_cppad
#
# initialize return variable
ok = True
# ---------------------------------------------------------------------
# number of dependent and independent variables
n = 3
#
# 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)
#
# ay[i] = j * ax[j] * ax[i]
# 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)
ax_j = ax[j]
ax_i = ax[i]
ay[i] = aj * ax_j * ax_i
#
#
# define af corresponding to f(x)
af = py_cppad.a_fun(ax, ay)
#
# Set select_d (domain) to all true,
# initial select_r (range) to all false
# initialize r to all zeros
select_d = py_cppad.vec_bool(n)
select_r = py_cppad.vec_bool(n)
r = py_cppad.vec_double(n)
for i in range( n ) :
select_d[i] = True
select_r[i] = False
r[i] = 0.0
#
#
# only select component n-1 of the range function
# f_0 (x) = (n-1) * x_{n-1} * x_0
select_r[0] = True
r[0] = 1.0
#
# sparisty pattern for Hessian
pattern = py_cppad.sparse_rc()
af.for_hes_sparsity(select_d, select_r, pattern)
#
# compute all possibly non-zero entries in Hessian
# (should only compute lower triangle becuase matrix is symmetric)
subset = py_cppad.sparse_rcv(pattern)
#
# work space used to save time for multiple calls
work = py_cppad.sparse_hes_work()
#
af.sparse_hes(subset, x, r, pattern, work)
#
# check that result is sparsity pattern for Hessian of f_0 (x)
ok = ok and subset.nnz() == 2
row = subset.row()
col = subset.col()
val = subset.val()
for k in range( 2 ) :
i = row[k]
j = col[k]
v = val[k]
if i <= j :
ok = ok and i == 0
ok = ok and j == n-1
#
if i >= j :
ok = ok and i == n-1
ok = ok and j == 0
#
ok = ok and v == n-1
#
#
return( ok )
#