Changeset 3811


Ignore:
Timestamp:
Mar 27, 2016 9:00:39 AM (4 years ago)
Author:
bradbell
Message:

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: f77cd2e443b8d14f354f97efc4c26369227dfe81
end hash code: 33d1f5b8d837d83c8baf656045d85d8b45f0b297

commit 33d1f5b8d837d83c8baf656045d85d8b45f0b297
Author: Brad Bell <bradbell@…>
Date: Sun Mar 27 04:01:14 2016 -0700

  1. Test non-symmetric reverse mode in eigen_mat_inv (and fix this example).
  2. Advance version to cppad-20160327.

commit 238e214ff36ca835efca615c609576dc8bf5038d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 18:10:59 2016 -0700

eigen_mat_inv.hpp: use inverse matrix (since we are computing it).
eigen_mat_mul.hpp: comment seperator.
eigen_mat_inv.cpp: matrix during recording is non-singular.
eigen_mat_inv.cpp: test first and second order forward.

commit fc918b0476cc8ea66abdf2904a71fc93472d279d
Author: Brad Bell <bradbell@…>
Date: Sat Mar 26 16:26:51 2016 -0700

Add test_more/eigen_mat_inv.cpp.


eigen_mat_inv.hpp: fix calculation of vy.
eigen_mat_inv.cpp: fix section title.
CMakeLists.txt: no longer necessary to have special library for eigen tests.
eigen_mat_inv.cpp: Test with non-symmetric matrix.

Location:
trunk
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed
  • trunk/AUTHORS

    r3810 r3811  
    22             ===========================================
    33
    4 To date, 2016-03-26, Bradley M. Bell is the sole author of CppAD.
     4To date, 2016-03-27, Bradley M. Bell is the sole author of CppAD.
    55While Bradley M. Bell worked for the University of Washington during
    66the development of CppAD, the following are also true:
  • trunk/CMakeLists.txt

    r3810 r3811  
    1818
    1919# cppad_version is used by set_version.sh to get the version number.
    20 SET(cppad_version      "20160326" )
     20SET(cppad_version      "20160327" )
    2121SET(cppad_url          "http://www.coin-or.org/CppAD" )
    2222SET(cppad_description  "Differentiation of C++ Algorithms" )
  • trunk/configure

    r3810 r3811  
    11#! /bin/sh
    22# Guess values for system-dependent variables and create Makefiles.
    3 # Generated by GNU Autoconf 2.69 for cppad 20160326.
     3# Generated by GNU Autoconf 2.69 for cppad 20160327.
    44#
    55# Report bugs to <cppad@list.coin-or.org>.
     
    581581PACKAGE_NAME='cppad'
    582582PACKAGE_TARNAME='cppad'
    583 PACKAGE_VERSION='20160326'
    584 PACKAGE_STRING='cppad 20160326'
     583PACKAGE_VERSION='20160327'
     584PACKAGE_STRING='cppad 20160327'
    585585PACKAGE_BUGREPORT='cppad@list.coin-or.org'
    586586PACKAGE_URL=''
     
    14071407  # This message is too long to be a string in the A/UX 3.1 sh.
    14081408  cat <<_ACEOF
    1409 \`configure' configures cppad 20160326 to adapt to many kinds of systems.
     1409\`configure' configures cppad 20160327 to adapt to many kinds of systems.
    14101410
    14111411Usage: $0 [OPTION]... [VAR=VALUE]...
     
    14771477if test -n "$ac_init_help"; then
    14781478  case $ac_init_help in
    1479      short | recursive ) echo "Configuration of cppad 20160326:";;
     1479     short | recursive ) echo "Configuration of cppad 20160327:";;
    14801480   esac
    14811481  cat <<\_ACEOF
     
    16111611if $ac_init_version; then
    16121612  cat <<\_ACEOF
    1613 cppad configure 20160326
     1613cppad configure 20160327
    16141614generated by GNU Autoconf 2.69
    16151615
     
    22402240running configure, to aid debugging if configure makes a mistake.
    22412241
    2242 It was created by cppad $as_me 20160326, which was
     2242It was created by cppad $as_me 20160327, which was
    22432243generated by GNU Autoconf 2.69.  Invocation command line was
    22442244
     
    31303130# Define the identity of the package.
    31313131 PACKAGE='cppad'
    3132  VERSION='20160326'
     3132 VERSION='20160327'
    31333133
    31343134
     
    85648564# values after options handling.
    85658565ac_log="
    8566 This file was extended by cppad $as_me 20160326, which was
     8566This file was extended by cppad $as_me 20160327, which was
    85678567generated by GNU Autoconf 2.69.  Invocation command line was
    85688568
     
    86218621ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
    86228622ac_cs_version="\\
    8623 cppad config.status 20160326
     8623cppad config.status 20160327
    86248624configured by $0, generated by GNU Autoconf 2.69,
    86258625  with options \\"\$ac_cs_config\\"
  • trunk/configure.ac

    r3810 r3811  
    1313dnl Process this file with autoconf to produce a configure script.
    1414dnl   package   version              bug-report
    15 AC_INIT([cppad], [20160326], [cppad@list.coin-or.org])
     15AC_INIT([cppad], [20160327], [cppad@list.coin-or.org])
    1616AM_SILENT_RULES([yes])
    1717
  • trunk/cppad/example/eigen_mat_inv.hpp

    r3810 r3811  
    3131$head Theory$$
    3232
    33 $head Forward$$
    34 For $latex k = 0 , \ldots$$, the $th k$$ order Taylor coefficient is given by
    35 $latex \[
    36         I_k = \sum_{\ell=0}^k A_\ell R_{k-\ell}
    37 \] $$
    38 $latex \[
    39         R_k = A_0^{-1} \left( I_k - \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
    40 \] $$
    41 where $latex I_k$$ is the $th k$$ order Taylor coefficient for the
    42 identity matrix. i.e., $latex I_k$$ is the identity matrix if $latex k = 0$$
    43 and is the zero matrix if $latex k \neq 0$$.
    44 
    45 $head Reverse$$
     33$subhead Forward$$
     34The zero order forward mode Taylor coefficient is give by
     35$latex \[
     36        R_0 = A_0^{-1}
     37\]$$
     38For $latex k = 1 , \ldots$$, the $th k$$ order Taylor coefficient is given by
     39$latex \[
     40        R_k = - R_0 \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
     41\] $$
     42
     43
     44$subhead Product of Three Matrices$$
     45Suppose $latex D = A B C$$, and $latex \bar{D}$$ is the partial of the
     46scalar final result with respect to $latex D$$. It follows that
     47$latex \[
     48        \R{d} D = \R{d} A B C + A \R{d} B C +  A B \R{d} C
     49\] $$
     50$latex \[
     51        \R{tr} ( \bar{D}^\R{T} \R{d} D )
     52        =
     53        \R{tr} ( \bar{D}^\R{T} \R{d} A B C ) +
     54        \R{tr} ( \bar{D}^\R{T} A \R{d} B C ) +
     55        \R{tr} ( \bar{D}^\R{T} A B \R{d} C )
     56\] $$
     57$latex \[
     58        \R{tr} ( \bar{D}^\R{T} \R{d} D )
     59        =
     60        \R{tr} ( B C \bar{D}^\R{T} \R{d} A ) +
     61        \R{tr} ( C \bar{D}^\R{T} A \R{d} B ) +
     62        \R{tr} ( \bar{D}^\R{T} A B \R{d} C )
     63\] $$
     64$latex \[
     65        \bar{A} = \bar{D} (B C)^\R{T} \W{,}
     66        \bar{B} = A^\R{T} \bar{D} C^\R{T} \W{,}
     67        \bar{C} = (A B)^\R{T} \bar{D}
     68\] $$
     69
     70$subhead Reverse$$
    4671We use $latex \bar{R}_k$$ for the partial of the scalar final result
    4772with respect to $latex R_k$$.
    48 Note that $latex R_0 = A_0^{-1}$$.
    4973The back-propagation algorithm that eliminates $latex R_k$$,
    5074for $latex k > 0$$, is
    51 for $latex \ell = 1 , \ldots , k$$
    52 $latex \[
    53 \bar{A}_\ell = \bar{A}_\ell - R_0 \bar{R}_k R_{k-\ell}^\R{T}
    54 \] $$
    55 $latex \[
    56 \bar{R}_{k-\ell} = \bar{R}_{k-\ell} - R_0 A_\ell^\R{T} \bar{R}_k
    57 \] $$
    58 $latex \[
    59 \bar{R}_0  = \bar{R}_0 + \bar{R}_k
    60         \left( I_k - \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T}
     75$latex \[
     76\bar{R}_0  = \bar{R}_0 - \bar{R}_k
     77        \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T}
    6178\] $$
    6279$latex \[
    6380\bar{R}_0  = \bar{R}_0 + \bar{R}_k ( A_0 R_k )^\R{T}
     81\] $$
     82and for $latex \ell = 1 , \ldots , k$$
     83$latex \[
     84\bar{A}_\ell = \bar{A}_\ell - R_0^\R{T} \bar{R}_k R_{k-\ell}^\R{T}
     85\] $$
     86$latex \[
     87\bar{R}_{k-\ell} = \bar{R}_{k-\ell} - ( R_0 A_\ell )^\R{T} \bar{R}_k
    6488\] $$
    6589The back-propagation algorithm that eliminates $latex R_0$$ is
     
    223247                // result for each order
    224248                //
    225                 // compute LU factorixation of arg_[0]
    226                 Eigen::FullPivLU<matrix> lu_arg( f_arg_[0] );
    227                 for(size_t k = 0; k < n_order; k++)
     249                f_result_[0] = f_arg_[0].inverse();
     250                for(size_t k = 1; k < n_order; k++)
    228251                {       // initialize sum
    229                         if( k == 0 )
    230                                 f_sum_ = matrix::Identity(nr_, nr_);
    231                         else
    232                                 f_sum_ = matrix::Zero(nr_, nr_);
     252                        f_sum_ = matrix::Zero(nr_, nr_);
    233253                        // compute sum
    234254                        for(size_t ell = 1; ell <= k; ell++)
    235255                                f_sum_ -= f_arg_[ell] * f_result_[k-ell];
    236                         // solve arg_[0] * result_[k] = sum
    237                         // result_[k] = arg_[0]^{-1} * sum
    238                         f_result_[k] = lu_arg.solve(f_sum_);
     256                        // result_[k] = arg_[0]^{-1} * sum_
     257                        f_result_[k] = f_result_[0] * f_sum_;
    239258                }
    240259                // -------------------------------------------------------------------
     
    244263                                ty[ i * n_order + k ] = f_result_[k].data()[i];
    245264                }
    246 
     265                // -------------------------------------------------------------------
    247266                // check if we are computing vy
    248267                if( vx.size() == 0 )
     
    254273                assert( n_order == 1 );
    255274                for(size_t j = 0; j < nr_; j++)
    256                 {       for(size_t i = 0; i < nr_; i++)
     275                {       // only row j of this column of the identity is non-zero
     276                        // initialize which elements of column j of result are variables
     277                        for(size_t i = 0; i < nr_; i++)
    257278                        {       // initialize vy as false
    258279                                size_t index = i * nr_ + j;
    259280                                vy[index]    = false;
    260281                        }
    261                         for(size_t i = 0; i < nr_; i++)
     282                        // determine if any elements in row j of argument are variables
     283                        bool row_var = false;
     284                        for(size_t ell = 0; ell < nr_; ell++)
     285                        {       // arg information
     286                                size_t index  = j * nr_ + ell;
     287                                row_var      |= vx[index];
     288                        }
     289                        if( row_var )
    262290                        {       for(size_t ell = 0; ell < nr_; ell++)
    263291                                {       // arg information
    264                                         size_t index   = i * nr_ + ell;
    265                                         bool var_arg  = vx[index];
    266                                         // identity information
    267                                         bool nz_eye   = ell == j;
    268                                         // result information
    269                                         index = i * nr_ + j;
    270                                         vy[index] |= var_arg & nz_eye;
     292                                        size_t index = j * nr_ + ell;
     293                                        bool not_zero = f_arg_[0](j, ell) != scalar(0.0);
     294                                        bool var      = vx[index];
     295                                        if( not_zero | var )
     296                                        {       // result information
     297                                                index = ell * nr_ + j;
     298                                                vy[index] = true;
     299                                        }
    271300                                }
    272301                        }
     
    336365                for(size_t k1 = n_order; k1 > 1; k1--)
    337366                {       size_t k = k1 - 1;
    338                         for(size_t ell = 1; ell <= k; ell++)
    339                         {       r_arg_[ell]    -=
    340                                         f_result_[0] * r_result_[k] * f_result_[k-ell].transpose();
    341                                 //
    342                                 r_result_[k-ell] -=
    343                                         f_result_[0] * f_arg_[ell].transpose() * r_result_[k];
    344                         }
     367                        // bar{R}_0 = bar{R}_0 + bar{R}_k (A_0 R_k)^T
    345368                        r_result_[0] +=
    346369                        r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose();
     370                        //
     371                        for(size_t ell = 1; ell <= k; ell++)
     372                        {       // bar{A}_l = bar{A}_l - R_0^T bar{R}_k R_{k-l}^T
     373                                r_arg_[ell] -= f_result_[0].transpose()
     374                                        * r_result_[k] * f_result_[k-ell].transpose();
     375                                // bar{R}_{k-l} = bar{R}_{k-1} - (R_0 A_l)^T bar{R}_k
     376                                r_result_[k-ell] -= f_arg_[ell].transpose()
     377                                        * f_result_[0].transpose() * r_result_[k];
     378                        }
    347379                }
    348380                r_arg_[0] -=
  • trunk/cppad/example/eigen_mat_mul.hpp

    r3810 r3811  
    250250                                ty[ i * n_order + k ] = f_result_[k].data()[i];
    251251                }
    252 
     252                // ------------------------------------------------------------------
    253253                // check if we are computing vy
    254254                if( vx.size() == 0 )
  • trunk/doc.omh

    r3810 r3811  
    9292$comment bin/version assumes that : follows cppad version number here$$
    9393$section
    94 cppad-20160326: A Package for Differentiation of C++ Algorithms
     94cppad-20160327: A Package for Differentiation of C++ Algorithms
    9595$$
    9696$mindex AD algorithmic differentiation automatic C++ algorithm derivative CppAD version cppad.hpp$$
  • trunk/example/atomic/eigen_mat_inv.cpp

    r3810 r3811  
    1818$$
    1919
    20 $section  Atomic Eigen Matrix Division: Example and Test$$
     20$section  Atomic Eigen Matrix Inverse: Example and Test$$
    2121
    2222$head Description$$
     
    8484        CPPAD_TESTVECTOR(ad_scalar) ad_x(n);
    8585        for(size_t j = 0; j < n; j++)
    86                 ad_x[j] = ad_scalar(j);
     86                ad_x[j] = ad_scalar(j + 1);
    8787        CppAD::Independent(ad_x);
    8888        // -------------------------------------------------------------------
  • trunk/omh/install/download.omh

    r3810 r3811  
    9797$rnext
    9898current  $cnext EPL $cnext $href%
    99 http://www.coin-or.org/download/source/CppAD/cppad-20160326.epl.tgz%
    100 cppad-20160326.epl.tgz%$$
     99http://www.coin-or.org/download/source/CppAD/cppad-20160327.epl.tgz%
     100cppad-20160327.epl.tgz%$$
    101101$rnext
    102102current  $cnext GPL $cnext $href%
    103 http://www.coin-or.org/download/source/CppAD/cppad-20160326.gpl.tgz%
    104 cppad-20160326.gpl.tgz%$$
     103http://www.coin-or.org/download/source/CppAD/cppad-20160327.gpl.tgz%
     104cppad-20160327.gpl.tgz%$$
    105105$tend
    106106
  • trunk/omh/whats_new/whats_new_16.omh

    r3810 r3811  
    2929        xam
    3030        makefile
     31        vy
    3132$$
    3233
     
    3940assist you in learning about changes between various versions of CppAD.
    4041
     42$head 03-27$$
     43Fix a bug in the calculation of the $cref atomic_eigen_mat_inv.hpp$$
     44$cref/reverse/atomic_eigen_mat_inv.hpp/Private/reverse/$$ example.
     45
    4146$head 03-26$$
     47$list number$$
    4248Implement and test the $cref atomic_eigen_mat_inv.cpp$$
    4349$cref/reverse/atomic_eigen_mat_inv.hpp/Private/reverse/$$ is implemented.
     50$lnext
     51Fix a bug in the calculation of
     52$cref/vy/atomic_forward/vy/$$ in the $cref atomic_eigen_mat_inv.hpp$$
     53$cref/forward/atomic_eigen_mat_inv.hpp/Private/forward/$$ example.
     54$lend
     55
    4456
    4557$head 03-25$$
  • trunk/test_more/CMakeLists.txt

    r3803 r3811  
    2121# ADD_DEFINITIONS("-DCPPAD_MAX_NUM_THREADS=1")
    2222
    23 # adolc_prefix
     23# adolc_sources, adolc_libs and CPPAD_ADOLC_TEST
    2424SET(sources base_adolc.cpp)
    2525sources_libs_define(adolc "${sources}" adolc ADOLC_TEST)
     26
     27# eigen_sources and CPPAD_EIGEN_TEST
     28SET(sources cppad_eigen.cpp eigen_mat_inv.cpp)
     29sources_libs_define(eigen "${sources}" "" EIGEN_TEST)
    2630
    2731# ipopt_prefix
     
    2933sources_libs_define(ipopt "${sources}" "${ipopt_LIBRARIES}" IPOPT_TEST)
    3034
    31 # sources that use eigen
    32 IF ( eigen_prefix )
    33         # compile eigen library separately so can use different compiler flags
    34         ADD_LIBRARY(test_more_eigen_lib EXCLUDE_FROM_ALL
    35                 cppad_eigen.cpp
    36         )
    37         # Adds -D define flags to the compilation of source files.
    38         ADD_DEFINITIONS("-DCPPAD_EIGEN_TEST")
    39         # Add other compiler flags
    40         add_cppad_cxx_flags(test_more_eigen_lib)
    41         #
    42         # Add eigen to list of libraries
    43         SET(eigen_libs test_more_eigen_lib)
    44 ELSE ( eigen_prefix )
    45         SET(eigen_sources "")
    46         SET(eigen_libs "")
    47 ENDIF ( eigen_prefix )
    48 
    49 # add_executable(<name> [WIN32] [MACOSX_BUNDLE] [EXCLUDE_FROM_ALL]
    50 #                 source1 source2 ... sourceN
    5135# )
    5236ADD_EXECUTABLE(test_more EXCLUDE_FROM_ALL test_more.cpp
     37        ${adolc_sources}
     38        ${eigen_sources}
     39        ${ipopt_sources}
    5340        old_usead_2.cpp
    5441        old_usead_1.cpp
     
    5643        old_reciprocal.cpp
    5744        old_mat_mul.cpp
    58         ${adolc_sources}
    59         ${ipopt_sources}
    6045        abs.cpp
    6146        acos.cpp
     
    161146        ${adolc_libs}
    162147        ${ipopt_libs}
    163         ${eigen_libs}
    164148        ${colpack_libs}
    165149)
  • trunk/test_more/makefile.am

    r3803 r3811  
    5151EIGEN_EXTRA_FILES   =
    5252noinst_LIBRARIES    = libeigen.a
    53 libeigen_a_SOURCES  =  cppad_eigen.cpp
     53libeigen_a_SOURCES  =  cppad_eigen.cpp eigen_mat_mul_inv.cpp
    5454EIGEN_LIB           = -L. -leigen
    5555libeigen_a_CXXFLAGS = \
     
    5959         -I$(EIGEN_DIR)/include
    6060else
    61 EIGEN_EXTRA_FILES   = cppad_eigen.cpp
     61EIGEN_EXTRA_FILES   = cppad_eigen.cpp eigen_mat_mul_inv.cpp
    6262EIGEN_LIB           =
    6363endif
  • trunk/test_more/test_more.cpp

    r3793 r3811  
    5050extern bool DivEq(void);
    5151extern bool DivZeroOne(void);
     52extern bool eigen_mat_inv(void);
    5253extern bool erf(void);
    5354extern bool Exp(void);
     
    257258# endif
    258259# ifdef CPPAD_EIGEN_TEST
    259         ok &= Run( cppad_eigen, "cppad_eigen" );
     260        ok &= Run( cppad_eigen,     "cppad_eigen"    );
     261        ok &= Run( eigen_mat_inv,   "eigen_mat_inv"  );
    260262# endif
    261263# if ! CPPAD_EIGENVECTOR
Note: See TracChangeset for help on using the changeset viewer.