Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library sparse jac_sparsity sparse_jac_pattern_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

jac_sparsity-> sparse_jac_pattern_xam.cpp sparse_jac_pattern_xam.m sparse_jac_pattern_xam.pm sparse_jac_pattern_xam.py

sparse_jac_pattern_xam.pm

Headings

package sparse_jac_pattern_xam;
sub sparse_jac_pattern_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);
     #
     # create dependent variables ay with ay[i] = 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 $ay_i = $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_in = new pm_cppad::sparse_rc();
     $pat_in->resize($n, $n, $n);
     for(my $k = 0; $k < $n; $k++) {
          $pat_in->put($k, $k, $k);
     }
     #
     # loop over forward and reverse mode
     for(my $mode = 0; $mode < 2; $mode++) {
          my $pat_out = new pm_cppad::sparse_rc();
          if( $mode == 0  ) {
               $af->for_jac_sparsity($pat_in, $pat_out);
          }
          if( $mode == 1  ) {
               $af->rev_jac_sparsity($pat_in, $pat_out);
          }
          #
          # check that result is sparsity pattern for Jacobian
          $ok = $ok && $n == $pat_out->nnz();
          my $col_major = $pat_out->col_major();
          my $row = $pat_out->row();
          my $col = $pat_out->col();
          for(my $k = 0; $k < $n; $k++) {
               my $ell = $col_major->get($k);
               my $r = $row->get($ell);
               my $c = $col->get($ell);
               my $i = $c+1;
               if( $i >=  $n  ) {
                    $i = $i - $n;
               }
               $ok = $ok && $c == $k;
               $ok = $ok && $r == $i;
          }
     }
     #
     return( $ok );
}