Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse sparse_jac sparse_jac_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_jac-> sparse_jac_xam.cpp sparse_jac_xam.m sparse_jac_xam.pm sparse_jac_xam.py

sparse_jac_xam.pm

Headings

package sparse_jac_xam;
sub sparse_jac_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;
     # one
     my $aone = new pm_cppad::a_double(1.0);
     #
     # 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);
     #
     # create dependent variables ay with ay[i] = (j+1) * ax[j]
     # 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 $ay_i = ($aj + $aone) * $ax->get($j);
          $ay->set($i, $ay_i);
     }
     #
     # define af corresponding to f(x)
     my $af = new pm_cppad::a_fun($ax, $ay);
     #
     # sparsity pattern for identity matrix
     my $pat_eye = new pm_cppad::sparse_rc();
     $pat_eye->resize($n, $n, $n);
     for(my $k = 0; $k < $n; $k++) {
          $pat_eye->put($k, $k, $k);
     }
     #
     # sparsity pattern for the Jacobian
     my $pat_jac = new pm_cppad::sparse_rc();
     $af->for_jac_sparsity($pat_eye, $pat_jac);
     #
     # loop over forward and reverse mode
     for(my $mode = 0; $mode < 2; $mode++) {
          # compute all possibly non-zero entries in Jacobian
          my $subset = new pm_cppad::sparse_rcv($pat_jac);
          # work space used to save time for multiple calls
          my $work = new pm_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 && $n == $subset->nnz();
          my $col_major = $subset->col_major();
          my $row = $subset->row();
          my $col = $subset->col();
          my $val = $subset->val();
          for(my $k = 0; $k < $n; $k++) {
               my $ell = $col_major->get($k);
               my $r = $row->get($ell);
               my $c = $col->get($ell);
               my $v = $val->get($ell);
               my $i = $c+1;
               if( $i >=  $n  ) {
                    $i = $i - $n;
               }
               $ok = $ok && $c == $k;
               $ok = $ok && $r == $i;
               $ok = $ok && $v == $c + 1.0;
          }
     }
     #
     return( $ok );
}