Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse sparse_hes sparse_hes_xam.pm

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_hes-> sparse_hes_xam.cpp sparse_hes_xam.m sparse_hes_xam.pm sparse_hes_xam.py

sparse_hes_xam.pm

Headings

package sparse_hes_xam;
sub sparse_hes_xam() {
     # check for standard perl programming conventions
     use strict;
     use warnings;
     #
     # load the Cppad Swig library
     use pm_cppad;
     #
     # initilaize return variable
     my $ok = 1;
     # ---------------------------------------------------------------------
     # number of dependent and independent variables
     my $n = 3;
     #
     # create the independent variables ax
     my $x = new pm_cppad::vec_double($n);
     for(my $i = 0; $i < $n ; $i++) {
          $x->set($i, $i + 2.0);
     }
     my $ax = pm_cppad::independent($x);
     #
     # ay[i] = j * ax[j] * ax[i]
     # where i = mod(j + 1, n)
     my $ay = new pm_cppad::vec_a_double($n);
     for(my $j = 0; $j < $n ; $j++) {
          my $i = $j+1;
          if( $i >= $n  ) {
               $i = $i - $n;
          }
          my $aj = new pm_cppad::a_double($j);
          my $ax_j = $ax->get($j);
          my $ax_i = $ax->get($i);
          $ay->set($i, $aj * $ax_j * $ax_i);
     }
     #
     # define af corresponding to f(x)
     my $af = new pm_cppad::a_fun($ax, $ay);
     #
     # Set select_d (domain) to all true,
     # initial select_r (range) to all false
     # initialize r to all zeros
     my $select_d = new pm_cppad::vec_bool($n);
     my $select_r = new pm_cppad::vec_bool($n);
     my $r = new pm_cppad::vec_double($n);
     for(my $i = 0; $i < $n; $i++) {
          $select_d->set($i, 1);
          $select_r->set($i, 0);
          $r->set($i, 0.0);
     }
     #
     # only select component n-1 of the range function
     # f_0 (x) = (n-1) * x_{n-1} * x_0
     $select_r->set(0, 1);
     $r->set(0, 1.0);
     #
     # sparisty pattern for Hessian
     my $pattern = new pm_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)
     my $subset = new pm_cppad::sparse_rcv($pattern);
     #
     # work space used to save time for multiple calls
     my $work = new pm_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 && $subset->nnz() == 2 ;
     my $row = $subset->row();
     my $col = $subset->col();
     my $val = $subset->val();
     for(my $k = 0; $k < 2; $k++) {
          my $i = $row->get($k);
          my $j = $col->get($k);
          my $v = $val->get($k);
          if( $i <= $j  ) {
               $ok = $ok && $i == 0;
               $ok = $ok && $j == $n-1;
          }
          if( $i >= $j  ) {
               $ok = $ok && $i == $n-1;
               $ok = $ok && $j == 0;
          }
          $ok = $ok && $v == $n-1;
     }
     #
     return( $ok );
}