mirror of
https://github.com/firestar5683/StarPilot.git
synced 2026-07-19 00:02:08 +08:00
external folder
This commit is contained in:
George Hotz
+364
@@ -0,0 +1,364 @@
|
||||
# ifndef CPPAD_EXAMPLE_BASE_ADOLC_HPP
|
||||
# define CPPAD_EXAMPLE_BASE_ADOLC_HPP
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
/*
|
||||
$begin base_adolc.hpp$$
|
||||
$spell
|
||||
stringstream
|
||||
struct
|
||||
string
|
||||
setprecision
|
||||
str
|
||||
valgrind
|
||||
azmul
|
||||
expm1
|
||||
atanh
|
||||
acosh
|
||||
asinh
|
||||
erf
|
||||
ifndef
|
||||
define
|
||||
endif
|
||||
Rel
|
||||
codassign
|
||||
eps
|
||||
std
|
||||
abs_geq
|
||||
fabs
|
||||
cppad.hpp
|
||||
undef
|
||||
Lt
|
||||
Le
|
||||
Eq
|
||||
Ge
|
||||
Gt
|
||||
namespace
|
||||
cassert
|
||||
condassign
|
||||
hpp
|
||||
bool
|
||||
const
|
||||
Adolc
|
||||
adouble
|
||||
CondExpOp
|
||||
inline
|
||||
enum
|
||||
CppAD
|
||||
pow
|
||||
acos
|
||||
asin
|
||||
atan
|
||||
cos
|
||||
cosh
|
||||
exp
|
||||
sqrt
|
||||
$$
|
||||
|
||||
|
||||
$section Enable use of AD<Base> where Base is Adolc's adouble Type$$
|
||||
|
||||
$head Syntax$$
|
||||
$codei%# include <cppad/example/base_adolc.hpp>
|
||||
%$$
|
||||
$children%
|
||||
example/general/mul_level_adolc.cpp
|
||||
%$$
|
||||
|
||||
$head Example$$
|
||||
The file $cref mul_level_adolc.cpp$$ contains an example use of
|
||||
Adolc's $code adouble$$ type for a CppAD $icode Base$$ type.
|
||||
It returns true if it succeeds and false otherwise.
|
||||
The file $cref mul_level_adolc_ode.cpp$$ contains a more realistic
|
||||
(and complex) example.
|
||||
|
||||
$head Include Files$$
|
||||
This file $code base_adolc.hpp$$ requires $code adouble$$ to be defined.
|
||||
In addition, it is included before $code <cppad/cppad.hpp>$$,
|
||||
but it needs to include parts of CppAD that are used by this file.
|
||||
This is done with the following include commands:
|
||||
$srccode%cpp% */
|
||||
# include <adolc/adolc.h>
|
||||
# include <cppad/base_require.hpp>
|
||||
/* %$$
|
||||
|
||||
$head CondExpOp$$
|
||||
The type $code adouble$$ supports a conditional assignment function
|
||||
with the syntax
|
||||
$codei%
|
||||
condassign(%a%, %b%, %c%, %d%)
|
||||
%$$
|
||||
which evaluates to
|
||||
$codei%
|
||||
%a% = (%b% > 0) ? %c% : %d%;
|
||||
%$$
|
||||
This enables one to include conditionals in the recording of
|
||||
$code adouble$$ operations and later evaluation for different
|
||||
values of the independent variables
|
||||
(in the same spirit as the CppAD $cref CondExp$$ function).
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline adouble CondExpOp(
|
||||
enum CppAD::CompareOp cop ,
|
||||
const adouble &left ,
|
||||
const adouble &right ,
|
||||
const adouble &trueCase ,
|
||||
const adouble &falseCase )
|
||||
{ adouble result;
|
||||
switch( cop )
|
||||
{
|
||||
case CompareLt: // left < right
|
||||
condassign(result, right - left, trueCase, falseCase);
|
||||
break;
|
||||
|
||||
case CompareLe: // left <= right
|
||||
condassign(result, left - right, falseCase, trueCase);
|
||||
break;
|
||||
|
||||
case CompareEq: // left == right
|
||||
condassign(result, left - right, falseCase, trueCase);
|
||||
condassign(result, right - left, falseCase, result);
|
||||
break;
|
||||
|
||||
case CompareGe: // left >= right
|
||||
condassign(result, right - left, falseCase, trueCase);
|
||||
break;
|
||||
|
||||
case CompareGt: // left > right
|
||||
condassign(result, left - right, trueCase, falseCase);
|
||||
break;
|
||||
default:
|
||||
CppAD::ErrorHandler::Call(
|
||||
true , __LINE__ , __FILE__ ,
|
||||
"CppAD::CondExp",
|
||||
"Error: for unknown reason."
|
||||
);
|
||||
result = trueCase;
|
||||
}
|
||||
return result;
|
||||
}
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head CondExpRel$$
|
||||
The $cref/CPPAD_COND_EXP_REL/base_cond_exp/CondExpRel/$$ macro invocation
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
CPPAD_COND_EXP_REL(adouble)
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head EqualOpSeq$$
|
||||
The Adolc user interface does not specify a way to determine if
|
||||
two $code adouble$$ variables correspond to the same operations sequence.
|
||||
Make $code EqualOpSeq$$ an error if it gets used:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline bool EqualOpSeq(const adouble &x, const adouble &y)
|
||||
{ CppAD::ErrorHandler::Call(
|
||||
true , __LINE__ , __FILE__ ,
|
||||
"CppAD::EqualOpSeq(x, y)",
|
||||
"Error: adouble does not support EqualOpSeq."
|
||||
);
|
||||
return false;
|
||||
}
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head Identical$$
|
||||
The Adolc user interface does not specify a way to determine if an
|
||||
$code adouble$$ depends on the independent variables.
|
||||
To be safe (but slow) return $code false$$ in all the cases below.
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline bool IdenticalPar(const adouble &x)
|
||||
{ return false; }
|
||||
inline bool IdenticalZero(const adouble &x)
|
||||
{ return false; }
|
||||
inline bool IdenticalOne(const adouble &x)
|
||||
{ return false; }
|
||||
inline bool IdenticalEqualPar(const adouble &x, const adouble &y)
|
||||
{ return false; }
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head Integer$$
|
||||
$srccode%cpp% */
|
||||
inline int Integer(const adouble &x)
|
||||
{ return static_cast<int>( x.getValue() ); }
|
||||
/* %$$
|
||||
|
||||
$head azmul$$
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
CPPAD_AZMUL( adouble )
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head Ordered$$
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline bool GreaterThanZero(const adouble &x)
|
||||
{ return (x > 0); }
|
||||
inline bool GreaterThanOrZero(const adouble &x)
|
||||
{ return (x >= 0); }
|
||||
inline bool LessThanZero(const adouble &x)
|
||||
{ return (x < 0); }
|
||||
inline bool LessThanOrZero(const adouble &x)
|
||||
{ return (x <= 0); }
|
||||
inline bool abs_geq(const adouble& x, const adouble& y)
|
||||
{ return fabs(x) >= fabs(y); }
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head Unary Standard Math$$
|
||||
The following $cref/required/base_require/$$ functions
|
||||
are defined by the Adolc package for the $code adouble$$ base case:
|
||||
$pre
|
||||
$$
|
||||
$code acos$$,
|
||||
$code asin$$,
|
||||
$code atan$$,
|
||||
$code cos$$,
|
||||
$code cosh$$,
|
||||
$code exp$$,
|
||||
$code fabs$$,
|
||||
$code log$$,
|
||||
$code sin$$,
|
||||
$code sinh$$,
|
||||
$code sqrt$$,
|
||||
$code tan$$.
|
||||
|
||||
$head erf, asinh, acosh, atanh, expm1, log1p$$
|
||||
If the
|
||||
$cref/erf, asinh, acosh, atanh, expm1, log1p
|
||||
/base_std_math
|
||||
/erf, asinh, acosh, atanh, expm1, log1p
|
||||
/$$,
|
||||
functions are supported by the compiler,
|
||||
they must also be supported by a $icode Base$$ type;
|
||||
The adolc package does not support these functions so make
|
||||
their use an error:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
# define CPPAD_BASE_ADOLC_NO_SUPPORT(fun) \
|
||||
inline adouble fun(const adouble& x) \
|
||||
{ CPPAD_ASSERT_KNOWN( \
|
||||
false, \
|
||||
#fun ": adolc does not support this function" \
|
||||
); \
|
||||
return 0.0; \
|
||||
}
|
||||
# if CPPAD_USE_CPLUSPLUS_2011
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(erf)
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(asinh)
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(acosh)
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(atanh)
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(expm1)
|
||||
CPPAD_BASE_ADOLC_NO_SUPPORT(log1p)
|
||||
# endif
|
||||
# undef CPPAD_BASE_ADOLC_NO_SUPPORT
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head sign$$
|
||||
This $cref/required/base_require/$$ function is defined using the
|
||||
$code codassign$$ function so that its $code adouble$$ operation sequence
|
||||
does not depend on the value of $icode x$$.
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline adouble sign(const adouble& x)
|
||||
{ adouble s_plus, s_minus, half(.5);
|
||||
// set s_plus to sign(x)/2, except for case x == 0, s_plus = -.5
|
||||
condassign(s_plus, +x, -half, +half);
|
||||
// set s_minus to -sign(x)/2, except for case x == 0, s_minus = -.5
|
||||
condassign(s_minus, -x, -half, +half);
|
||||
// set s to sign(x)
|
||||
return s_plus - s_minus;
|
||||
}
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head abs$$
|
||||
This $cref/required/base_require/$$ function uses the adolc $code fabs$$
|
||||
function:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
inline adouble abs(const adouble& x)
|
||||
{ return fabs(x); }
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head pow$$
|
||||
This $cref/required/base_require/$$ function
|
||||
is defined by the Adolc package for the $code adouble$$ base case.
|
||||
|
||||
$head numeric_limits$$
|
||||
The following defines the CppAD $cref numeric_limits$$
|
||||
for the type $code adouble$$:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
CPPAD_NUMERIC_LIMITS(double, adouble)
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head to_string$$
|
||||
The following defines the CppAD $cref to_string$$ function
|
||||
for the type $code adouble$$:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
template <> struct to_string_struct<adouble>
|
||||
{ std::string operator()(const adouble& x)
|
||||
{ std::stringstream os;
|
||||
int n_digits = 1 + std::numeric_limits<double>::digits10;
|
||||
os << std::setprecision(n_digits);
|
||||
os << x.value();
|
||||
return os.str();
|
||||
}
|
||||
};
|
||||
}
|
||||
/* %$$
|
||||
|
||||
$head hash_code$$
|
||||
It appears that an $code adouble$$ object can have fields
|
||||
that are not initialized.
|
||||
This results in a $code valgrind$$ error when these fields are used by the
|
||||
$cref/default/base_hash/Default/$$ hashing function.
|
||||
For this reason, the $code adouble$$ class overrides the default definition.
|
||||
$srccode|cpp| */
|
||||
namespace CppAD {
|
||||
inline unsigned short hash_code(const adouble& x)
|
||||
{ unsigned short code = 0;
|
||||
double value = x.value();
|
||||
if( value == 0.0 )
|
||||
return code;
|
||||
double log_x = std::log( fabs( value ) );
|
||||
// assume log( std::numeric_limits<double>::max() ) is near 700
|
||||
code = static_cast<unsigned short>(
|
||||
(CPPAD_HASH_TABLE_SIZE / 700 + 1) * log_x
|
||||
);
|
||||
code = code % CPPAD_HASH_TABLE_SIZE;
|
||||
return code;
|
||||
}
|
||||
}
|
||||
/* |$$
|
||||
Note that after the hash codes match, the
|
||||
$cref/Identical/base_adolc.hpp/Identical/$$ function will be used
|
||||
to make sure two values are the same and one can replace the other.
|
||||
A more sophisticated implementation of the $code Identical$$ function
|
||||
would detect which $code adouble$$ values depend on the
|
||||
$code adouble$$ independent variables (and hence can change).
|
||||
|
||||
|
||||
$end
|
||||
*/
|
||||
# endif
|
||||
|
||||
@@ -0,0 +1,147 @@
|
||||
# ifndef CPPAD_EXAMPLE_CPPAD_EIGEN_HPP
|
||||
# define CPPAD_EXAMPLE_CPPAD_EIGEN_HPP
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
/*
|
||||
$begin cppad_eigen.hpp$$
|
||||
$spell
|
||||
impl
|
||||
typename
|
||||
Real Real
|
||||
inline
|
||||
neg
|
||||
eps
|
||||
atan
|
||||
Num
|
||||
acos
|
||||
asin
|
||||
CppAD
|
||||
std::numeric
|
||||
enum
|
||||
Mul
|
||||
Eigen
|
||||
cppad.hpp
|
||||
namespace
|
||||
struct
|
||||
typedef
|
||||
const
|
||||
imag
|
||||
sqrt
|
||||
exp
|
||||
cos
|
||||
$$
|
||||
$section Enable Use of Eigen Linear Algebra Package with CppAD$$
|
||||
|
||||
$head Syntax$$
|
||||
$codei%# include <cppad/example/cppad_eigen.hpp>%$$
|
||||
$children%
|
||||
example/general/eigen_array.cpp%
|
||||
example/general/eigen_det.cpp
|
||||
%$$
|
||||
|
||||
$head Purpose$$
|
||||
Enables the use of the $cref/eigen/eigen_prefix/$$
|
||||
linear algebra package with the type $icode%AD<%Base%>%$$; see
|
||||
$href%
|
||||
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html%
|
||||
custom scalar types
|
||||
%$$.
|
||||
|
||||
$head Example$$
|
||||
The files $cref eigen_array.cpp$$ and $cref eigen_det.cpp$$
|
||||
contain an example and test of this include file.
|
||||
They return true if they succeed and false otherwise.
|
||||
|
||||
$head Include Files$$
|
||||
The file $code cppad_eigen.hpp$$ includes both
|
||||
$code <cppad/cppad.hpp>$$ and $code <Eigen/Core>$$.
|
||||
$srccode%cpp% */
|
||||
# include <Eigen/Core>
|
||||
# include <cppad/cppad.hpp>
|
||||
/* %$$
|
||||
$head Eigen NumTraits$$
|
||||
Eigen needs the following definitions to work properly
|
||||
with $codei%AD<%Base%>%$$ scalars:
|
||||
$srccode%cpp% */
|
||||
namespace Eigen {
|
||||
template <class Base> struct NumTraits< CppAD::AD<Base> >
|
||||
{ // type that corresponds to the real part of an AD<Base> value
|
||||
typedef CppAD::AD<Base> Real;
|
||||
// type for AD<Base> operations that result in non-integer values
|
||||
typedef CppAD::AD<Base> NonInteger;
|
||||
// type to use for numeric literals such as "2" or "0.5".
|
||||
typedef CppAD::AD<Base> Literal;
|
||||
// type for nested value inside an AD<Base> expression tree
|
||||
typedef CppAD::AD<Base> Nested;
|
||||
|
||||
enum {
|
||||
// does not support complex Base types
|
||||
IsComplex = 0 ,
|
||||
// does not support integer Base types
|
||||
IsInteger = 0 ,
|
||||
// only support signed Base types
|
||||
IsSigned = 1 ,
|
||||
// must initialize an AD<Base> object
|
||||
RequireInitialization = 1 ,
|
||||
// computational cost of the corresponding operations
|
||||
ReadCost = 1 ,
|
||||
AddCost = 2 ,
|
||||
MulCost = 2
|
||||
};
|
||||
|
||||
// machine epsilon with type of real part of x
|
||||
// (use assumption that Base is not complex)
|
||||
static CppAD::AD<Base> epsilon(void)
|
||||
{ return CppAD::numeric_limits< CppAD::AD<Base> >::epsilon(); }
|
||||
|
||||
// relaxed version of machine epsilon for comparison of different
|
||||
// operations that should result in the same value
|
||||
static CppAD::AD<Base> dummy_precision(void)
|
||||
{ return 100. *
|
||||
CppAD::numeric_limits< CppAD::AD<Base> >::epsilon();
|
||||
}
|
||||
|
||||
// minimum normalized positive value
|
||||
static CppAD::AD<Base> lowest(void)
|
||||
{ return CppAD::numeric_limits< CppAD::AD<Base> >::min(); }
|
||||
|
||||
// maximum finite value
|
||||
static CppAD::AD<Base> highest(void)
|
||||
{ return CppAD::numeric_limits< CppAD::AD<Base> >::max(); }
|
||||
|
||||
// number of decimal digits that can be represented without change.
|
||||
static int digits10(void)
|
||||
{ return CppAD::numeric_limits< CppAD::AD<Base> >::digits10; }
|
||||
};
|
||||
}
|
||||
/* %$$
|
||||
$head CppAD Namespace$$
|
||||
Eigen also needs the following definitions to work properly
|
||||
with $codei%AD<%Base%>%$$ scalars:
|
||||
$srccode%cpp% */
|
||||
namespace CppAD {
|
||||
// functions that return references
|
||||
template <class Base> const AD<Base>& conj(const AD<Base>& x)
|
||||
{ return x; }
|
||||
template <class Base> const AD<Base>& real(const AD<Base>& x)
|
||||
{ return x; }
|
||||
|
||||
// functions that return values (note abs is defined by cppad.hpp)
|
||||
template <class Base> AD<Base> imag(const AD<Base>& x)
|
||||
{ return CppAD::AD<Base>(0.); }
|
||||
template <class Base> AD<Base> abs2(const AD<Base>& x)
|
||||
{ return x * x; }
|
||||
}
|
||||
|
||||
/* %$$
|
||||
$end
|
||||
*/
|
||||
# endif
|
||||
@@ -0,0 +1,377 @@
|
||||
// $Id$
|
||||
# ifndef CPPAD_EXAMPLE_EIGEN_CHOLESKY_HPP
|
||||
# define CPPAD_EXAMPLE_EIGEN_CHOLESKY_HPP
|
||||
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
|
||||
/*
|
||||
$begin atomic_eigen_cholesky.hpp$$
|
||||
$spell
|
||||
Eigen
|
||||
Taylor
|
||||
Cholesky
|
||||
op
|
||||
$$
|
||||
|
||||
$section Atomic Eigen Cholesky Factorization Class$$
|
||||
|
||||
$head Purpose$$
|
||||
Construct an atomic operation that computes a lower triangular matrix
|
||||
$latex L $$ such that $latex L L^\R{T} = A$$
|
||||
for any positive integer $latex p$$
|
||||
and symmetric positive definite matrix $latex A \in \B{R}^{p \times p}$$.
|
||||
|
||||
$head Start Class Definition$$
|
||||
$srccode%cpp% */
|
||||
# include <cppad/cppad.hpp>
|
||||
# include <Eigen/Dense>
|
||||
|
||||
|
||||
/* %$$
|
||||
$head Public$$
|
||||
|
||||
$subhead Types$$
|
||||
$srccode%cpp% */
|
||||
namespace { // BEGIN_EMPTY_NAMESPACE
|
||||
|
||||
template <class Base>
|
||||
class atomic_eigen_cholesky : public CppAD::atomic_base<Base> {
|
||||
public:
|
||||
// -----------------------------------------------------------
|
||||
// type of elements during calculation of derivatives
|
||||
typedef Base scalar;
|
||||
// type of elements during taping
|
||||
typedef CppAD::AD<scalar> ad_scalar;
|
||||
//
|
||||
// type of matrix during calculation of derivatives
|
||||
typedef Eigen::Matrix<
|
||||
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
|
||||
// type of matrix during taping
|
||||
typedef Eigen::Matrix<
|
||||
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
|
||||
//
|
||||
// lower triangular scalar matrix
|
||||
typedef Eigen::TriangularView<matrix, Eigen::Lower> lower_view;
|
||||
/* %$$
|
||||
$subhead Constructor$$
|
||||
$srccode%cpp% */
|
||||
// constructor
|
||||
atomic_eigen_cholesky(void) : CppAD::atomic_base<Base>(
|
||||
"atom_eigen_cholesky" ,
|
||||
CppAD::atomic_base<Base>::set_sparsity_enum
|
||||
)
|
||||
{ }
|
||||
/* %$$
|
||||
$subhead op$$
|
||||
$srccode%cpp% */
|
||||
// use atomic operation to invert an AD matrix
|
||||
ad_matrix op(const ad_matrix& arg)
|
||||
{ size_t nr = size_t( arg.rows() );
|
||||
size_t ny = ( (nr + 1 ) * nr ) / 2;
|
||||
size_t nx = 1 + ny;
|
||||
assert( nr == size_t( arg.cols() ) );
|
||||
// -------------------------------------------------------------------
|
||||
// packed version of arg
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
|
||||
size_t index = 0;
|
||||
packed_arg[index++] = ad_scalar( nr );
|
||||
// lower triangle of symmetric matrix A
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
packed_arg[index++] = arg(i, j);
|
||||
}
|
||||
assert( index == nx );
|
||||
// -------------------------------------------------------------------
|
||||
// packed version of result = arg^{-1}.
|
||||
// This is an atomic_base function call that CppAD uses to
|
||||
// store the atomic operation on the tape.
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
|
||||
(*this)(packed_arg, packed_result);
|
||||
// -------------------------------------------------------------------
|
||||
// unpack result matrix L
|
||||
ad_matrix result = ad_matrix::Zero(nr, nr);
|
||||
index = 0;
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
result(i, j) = packed_result[index++];
|
||||
}
|
||||
return result;
|
||||
}
|
||||
/* %$$
|
||||
$head Private$$
|
||||
|
||||
$subhead Variables$$
|
||||
$srccode%cpp% */
|
||||
private:
|
||||
// -------------------------------------------------------------
|
||||
// one forward mode vector of matrices for argument and result
|
||||
CppAD::vector<matrix> f_arg_, f_result_;
|
||||
// one reverse mode vector of matrices for argument and result
|
||||
CppAD::vector<matrix> r_arg_, r_result_;
|
||||
// -------------------------------------------------------------
|
||||
/* %$$
|
||||
$subhead forward$$
|
||||
$srccode%cpp% */
|
||||
// forward mode routine called by CppAD
|
||||
virtual bool forward(
|
||||
// lowest order Taylor coefficient we are evaluating
|
||||
size_t p ,
|
||||
// highest order Taylor coefficient we are evaluating
|
||||
size_t q ,
|
||||
// which components of x are variables
|
||||
const CppAD::vector<bool>& vx ,
|
||||
// which components of y are variables
|
||||
CppAD::vector<bool>& vy ,
|
||||
// tx [ j * (q+1) + k ] is x_j^k
|
||||
const CppAD::vector<scalar>& tx ,
|
||||
// ty [ i * (q+1) + k ] is y_i^k
|
||||
CppAD::vector<scalar>& ty
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
size_t ny = ((nr + 1) * nr) / 2;
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 1 + ny;
|
||||
# endif
|
||||
assert( vx.size() == 0 || nx == vx.size() );
|
||||
assert( vx.size() == 0 || ny == vy.size() );
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
//
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_arg_ and f_result_ are large enough
|
||||
assert( f_arg_.size() == f_result_.size() );
|
||||
if( f_arg_.size() < n_order )
|
||||
{ f_arg_.resize(n_order);
|
||||
f_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ f_arg_[k].resize(nr, nr);
|
||||
f_result_[k].resize(nr, nr);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_arg_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ size_t index = 1;
|
||||
// unpack arg values for this order
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
{ f_arg_[k](i, j) = tx[ index * n_order + k ];
|
||||
f_arg_[k](j, i) = f_arg_[k](i, j);
|
||||
index++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// result for each order
|
||||
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
|
||||
//
|
||||
Eigen::LLT<matrix> cholesky(f_arg_[0]);
|
||||
f_result_[0] = cholesky.matrixL();
|
||||
lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
|
||||
for(size_t k = 1; k < n_order; k++)
|
||||
{ // initialize sum as A_k
|
||||
matrix f_sum = f_arg_[k];
|
||||
// compute A_k - B_k
|
||||
for(size_t ell = 1; ell < k; ell++)
|
||||
f_sum -= f_result_[ell] * f_result_[k-ell].transpose();
|
||||
// compute L_0^{-1} * (A_k - B_k) * L_0^{-T}
|
||||
matrix temp = L_0.template solve<Eigen::OnTheLeft>(f_sum);
|
||||
temp = L_0.transpose().template solve<Eigen::OnTheRight>(temp);
|
||||
// divide the diagonal by 2
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
temp(i, i) /= scalar(2.0);
|
||||
// L_k = L_0 * low[ L_0^{-1} * (A_k - B_k) * L_0^{-T} ]
|
||||
lower_view view = temp.template triangularView<Eigen::Lower>();
|
||||
f_result_[k] = f_result_[0] * view;
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// pack result_ into ty
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ size_t index = 0;
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
{ ty[ index * n_order + k ] = f_result_[k](i, j);
|
||||
index++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// check if we are computing vy
|
||||
if( vx.size() == 0 )
|
||||
return true;
|
||||
// ------------------------------------------------------------------
|
||||
// This is a very dumb algorithm that over estimates which
|
||||
// elements of the inverse are variables (which is not efficient).
|
||||
bool var = false;
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
var |= vx[1 + i];
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
vy[i] = var;
|
||||
//
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead reverse$$
|
||||
$srccode%cpp% */
|
||||
// reverse mode routine called by CppAD
|
||||
virtual bool reverse(
|
||||
// highest order Taylor coefficient that we are computing derivative of
|
||||
size_t q ,
|
||||
// forward mode Taylor coefficients for x variables
|
||||
const CppAD::vector<double>& tx ,
|
||||
// forward mode Taylor coefficients for y variables
|
||||
const CppAD::vector<double>& ty ,
|
||||
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
|
||||
CppAD::vector<double>& px ,
|
||||
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
|
||||
const CppAD::vector<double>& py
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t ny = ( (nr + 1 ) * nr ) / 2;
|
||||
size_t nx = 1 + ny;
|
||||
# endif
|
||||
//
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
assert( px.size() == tx.size() );
|
||||
assert( py.size() == ty.size() );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_arg_ is large enough
|
||||
assert( f_arg_.size() == f_result_.size() );
|
||||
// must have previous run forward with order >= n_order
|
||||
assert( f_arg_.size() >= n_order );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure r_arg_, r_result_ are large enough
|
||||
assert( r_arg_.size() == r_result_.size() );
|
||||
if( r_arg_.size() < n_order )
|
||||
{ r_arg_.resize(n_order);
|
||||
r_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ r_arg_[k].resize(nr, nr);
|
||||
r_result_[k].resize(nr, nr);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_arg_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ size_t index = 1;
|
||||
// unpack arg values for this order
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
{ f_arg_[k](i, j) = tx[ index * n_order + k ];
|
||||
f_arg_[k](j, i) = f_arg_[k](i, j);
|
||||
index++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack py into r_result_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ r_result_[k] = matrix::Zero(nr, nr);
|
||||
size_t index = 0;
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j <= i; j++)
|
||||
{ r_result_[k](i, j) = py[ index * n_order + k ];
|
||||
index++;
|
||||
}
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// initialize r_arg_ as zero
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
r_arg_[k] = matrix::Zero(nr, nr);
|
||||
// -------------------------------------------------------------------
|
||||
// matrix reverse mode calculation
|
||||
lower_view L_0 = f_result_[0].template triangularView<Eigen::Lower>();
|
||||
//
|
||||
for(size_t k1 = n_order; k1 > 1; k1--)
|
||||
{ size_t k = k1 - 1;
|
||||
//
|
||||
// L_0^T * bar{L}_k
|
||||
matrix tmp1 = L_0.transpose() * r_result_[k];
|
||||
//
|
||||
//low[ L_0^T * bar{L}_k ]
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
tmp1(i, i) /= scalar(2.0);
|
||||
matrix tmp2 = tmp1.template triangularView<Eigen::Lower>();
|
||||
//
|
||||
// L_0^{-T} low[ L_0^T * bar{L}_k ]
|
||||
tmp1 = L_0.transpose().template solve<Eigen::OnTheLeft>( tmp2 );
|
||||
//
|
||||
// M_k = L_0^{-T} * low[ L_0^T * bar{L}_k ]^{T} L_0^{-1}
|
||||
matrix M_k = L_0.transpose().template
|
||||
solve<Eigen::OnTheLeft>( tmp1.transpose() );
|
||||
//
|
||||
// remove L_k and compute bar{B}_k
|
||||
matrix barB_k = scalar(0.5) * ( M_k + M_k.transpose() );
|
||||
r_arg_[k] += barB_k;
|
||||
barB_k = scalar(-1.0) * barB_k;
|
||||
//
|
||||
// 2.0 * lower( bar{B}_k L_k )
|
||||
matrix temp = scalar(2.0) * barB_k * f_result_[k];
|
||||
temp = temp.template triangularView<Eigen::Lower>();
|
||||
//
|
||||
// remove C_k
|
||||
r_result_[0] += temp;
|
||||
//
|
||||
// remove B_k
|
||||
for(size_t ell = 1; ell < k; ell++)
|
||||
{ // bar{L}_ell = 2 * lower( \bar{B}_k * L_{k-ell} )
|
||||
temp = scalar(2.0) * barB_k * f_result_[k-ell];
|
||||
r_result_[ell] += temp.template triangularView<Eigen::Lower>();
|
||||
}
|
||||
}
|
||||
// M_0 = L_0^{-T} * low[ L_0^T * bar{L}_0 ]^{T} L_0^{-1}
|
||||
matrix M_0 = L_0.transpose() * r_result_[0];
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
M_0(i, i) /= scalar(2.0);
|
||||
M_0 = M_0.template triangularView<Eigen::Lower>();
|
||||
M_0 = L_0.template solve<Eigen::OnTheRight>( M_0 );
|
||||
M_0 = L_0.transpose().template solve<Eigen::OnTheLeft>( M_0 );
|
||||
// remove L_0
|
||||
r_arg_[0] += scalar(0.5) * ( M_0 + M_0.transpose() );
|
||||
// -------------------------------------------------------------------
|
||||
// pack r_arg into px
|
||||
// note that only the lower triangle of barA_k is stored in px
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ size_t index = 0;
|
||||
px[ index * n_order + k ] = 0.0;
|
||||
index++;
|
||||
for(size_t i = 0; i < nr; i++)
|
||||
{ for(size_t j = 0; j < i; j++)
|
||||
{ px[ index * n_order + k ] = 2.0 * r_arg_[k](i, j);
|
||||
index++;
|
||||
}
|
||||
px[ index * n_order + k] = r_arg_[k](i, i);
|
||||
index++;
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head End Class Definition$$
|
||||
$srccode%cpp% */
|
||||
}; // End of atomic_eigen_cholesky class
|
||||
|
||||
} // END_EMPTY_NAMESPACE
|
||||
/* %$$
|
||||
$end
|
||||
*/
|
||||
|
||||
|
||||
# endif
|
||||
@@ -0,0 +1,396 @@
|
||||
// $Id$
|
||||
# ifndef CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
|
||||
# define CPPAD_EXAMPLE_EIGEN_MAT_INV_HPP
|
||||
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
|
||||
/*
|
||||
$begin atomic_eigen_mat_inv.hpp$$
|
||||
$spell
|
||||
Eigen
|
||||
Taylor
|
||||
$$
|
||||
|
||||
$section Atomic Eigen Matrix Inversion Class$$
|
||||
|
||||
$head Purpose$$
|
||||
Construct an atomic operation that computes the matrix inverse
|
||||
$latex R = A^{-1}$$
|
||||
for any positive integer $latex p$$
|
||||
and invertible matrix $latex A \in \B{R}^{p \times p}$$.
|
||||
|
||||
$head Matrix Dimensions$$
|
||||
This example puts the matrix dimension $latex p$$
|
||||
in the atomic function arguments,
|
||||
instead of the $cref/constructor/atomic_ctor/$$,
|
||||
so it can be different for different calls to the atomic function.
|
||||
|
||||
$head Theory$$
|
||||
|
||||
$subhead Forward$$
|
||||
The zero order forward mode Taylor coefficient is give by
|
||||
$latex \[
|
||||
R_0 = A_0^{-1}
|
||||
\]$$
|
||||
For $latex k = 1 , \ldots$$,
|
||||
the $th k$$ order Taylor coefficient of $latex A R$$ is given by
|
||||
$latex \[
|
||||
0 = \sum_{\ell=0}^k A_\ell R_{k-\ell}
|
||||
\] $$
|
||||
Solving for $latex R_k$$ in terms of the coefficients
|
||||
for $latex A$$ and the lower order coefficients for $latex R$$ we have
|
||||
$latex \[
|
||||
R_k = - R_0 \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
|
||||
\] $$
|
||||
Furthermore, once we have $latex R_k$$ we can compute the sum using
|
||||
$latex \[
|
||||
A_0 R_k = - \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)
|
||||
\] $$
|
||||
|
||||
|
||||
$subhead Product of Three Matrices$$
|
||||
Suppose $latex \bar{E}$$ is the derivative of the
|
||||
scalar value function $latex s(E)$$ with respect to $latex E$$; i.e.,
|
||||
$latex \[
|
||||
\bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} }
|
||||
\] $$
|
||||
Also suppose that $latex t$$ is a scalar valued argument and
|
||||
$latex \[
|
||||
E(t) = B(t) C(t) D(t)
|
||||
\] $$
|
||||
It follows that
|
||||
$latex \[
|
||||
E'(t) = B'(t) C(t) D(t) + B(t) C'(t) D(t) + B(t) C(t) D'(t)
|
||||
\] $$
|
||||
|
||||
$latex \[
|
||||
(s \circ E)'(t)
|
||||
=
|
||||
\R{tr} [ \bar{E}^\R{T} E'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
=
|
||||
\R{tr} [ \bar{E}^\R{T} B'(t) C(t) D(t) ] +
|
||||
\R{tr} [ \bar{E}^\R{T} B(t) C'(t) D(t) ] +
|
||||
\R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
=
|
||||
\R{tr} [ B(t) D(t) \bar{E}^\R{T} B'(t) ] +
|
||||
\R{tr} [ D(t) \bar{E}^\R{T} B(t) C'(t) ] +
|
||||
\R{tr} [ \bar{E}^\R{T} B(t) C(t) D'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
\bar{B} = \bar{E} (C D)^\R{T} \W{,}
|
||||
\bar{C} = B^\R{T} \bar{E} D^\R{T} \W{,}
|
||||
\bar{D} = (B C)^\R{T} \bar{E}
|
||||
\] $$
|
||||
|
||||
$subhead Reverse$$
|
||||
For $latex k > 0$$, reverse mode
|
||||
eliminates $latex R_k$$ and expresses the function values
|
||||
$latex s$$ in terms of the coefficients of $latex A$$
|
||||
and the lower order coefficients of $latex R$$.
|
||||
The effect on $latex \bar{R}_0$$
|
||||
(of eliminating $latex R_k$$) is
|
||||
$latex \[
|
||||
\bar{R}_0
|
||||
= \bar{R}_0 - \bar{R}_k \left( \sum_{\ell=1}^k A_\ell R_{k-\ell} \right)^\R{T}
|
||||
= \bar{R}_0 + \bar{R}_k ( A_0 R_k )^\R{T}
|
||||
\] $$
|
||||
For $latex \ell = 1 , \ldots , k$$,
|
||||
the effect on $latex \bar{R}_{k-\ell}$$ and $latex A_\ell$$
|
||||
(of eliminating $latex R_k$$) is
|
||||
$latex \[
|
||||
\bar{A}_\ell = \bar{A}_\ell - R_0^\R{T} \bar{R}_k R_{k-\ell}^\R{T}
|
||||
\] $$
|
||||
$latex \[
|
||||
\bar{R}_{k-\ell} = \bar{R}_{k-\ell} - ( R_0 A_\ell )^\R{T} \bar{R}_k
|
||||
\] $$
|
||||
We note that
|
||||
$latex \[
|
||||
R_0 '(t) A_0 (t) + R_0 (t) A_0 '(t) = 0
|
||||
\] $$
|
||||
$latex \[
|
||||
R_0 '(t) = - R_0 (t) A_0 '(t) R_0 (t)
|
||||
\] $$
|
||||
The reverse mode formula that eliminates $latex R_0$$ is
|
||||
$latex \[
|
||||
\bar{A}_0
|
||||
= \bar{A}_0 - R_0^\R{T} \bar{R}_0 R_0^\R{T}
|
||||
\]$$
|
||||
|
||||
$nospell
|
||||
|
||||
$head Start Class Definition$$
|
||||
$srccode%cpp% */
|
||||
# include <cppad/cppad.hpp>
|
||||
# include <Eigen/Core>
|
||||
# include <Eigen/LU>
|
||||
|
||||
|
||||
|
||||
/* %$$
|
||||
$head Public$$
|
||||
|
||||
$subhead Types$$
|
||||
$srccode%cpp% */
|
||||
namespace { // BEGIN_EMPTY_NAMESPACE
|
||||
|
||||
template <class Base>
|
||||
class atomic_eigen_mat_inv : public CppAD::atomic_base<Base> {
|
||||
public:
|
||||
// -----------------------------------------------------------
|
||||
// type of elements during calculation of derivatives
|
||||
typedef Base scalar;
|
||||
// type of elements during taping
|
||||
typedef CppAD::AD<scalar> ad_scalar;
|
||||
// type of matrix during calculation of derivatives
|
||||
typedef Eigen::Matrix<
|
||||
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
|
||||
// type of matrix during taping
|
||||
typedef Eigen::Matrix<
|
||||
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
|
||||
/* %$$
|
||||
$subhead Constructor$$
|
||||
$srccode%cpp% */
|
||||
// constructor
|
||||
atomic_eigen_mat_inv(void) : CppAD::atomic_base<Base>(
|
||||
"atom_eigen_mat_inv" ,
|
||||
CppAD::atomic_base<Base>::set_sparsity_enum
|
||||
)
|
||||
{ }
|
||||
/* %$$
|
||||
$subhead op$$
|
||||
$srccode%cpp% */
|
||||
// use atomic operation to invert an AD matrix
|
||||
ad_matrix op(const ad_matrix& arg)
|
||||
{ size_t nr = size_t( arg.rows() );
|
||||
size_t ny = nr * nr;
|
||||
size_t nx = 1 + ny;
|
||||
assert( nr == size_t( arg.cols() ) );
|
||||
// -------------------------------------------------------------------
|
||||
// packed version of arg
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
|
||||
packed_arg[0] = ad_scalar( nr );
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
packed_arg[1 + i] = arg.data()[i];
|
||||
// -------------------------------------------------------------------
|
||||
// packed version of result = arg^{-1}.
|
||||
// This is an atomic_base function call that CppAD uses to
|
||||
// store the atomic operation on the tape.
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
|
||||
(*this)(packed_arg, packed_result);
|
||||
// -------------------------------------------------------------------
|
||||
// unpack result matrix
|
||||
ad_matrix result(nr, nr);
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
result.data()[i] = packed_result[i];
|
||||
return result;
|
||||
}
|
||||
/* %$$
|
||||
$head Private$$
|
||||
|
||||
$subhead Variables$$
|
||||
$srccode%cpp% */
|
||||
private:
|
||||
// -------------------------------------------------------------
|
||||
// one forward mode vector of matrices for argument and result
|
||||
CppAD::vector<matrix> f_arg_, f_result_;
|
||||
// one reverse mode vector of matrices for argument and result
|
||||
CppAD::vector<matrix> r_arg_, r_result_;
|
||||
// -------------------------------------------------------------
|
||||
/* %$$
|
||||
$subhead forward$$
|
||||
$srccode%cpp% */
|
||||
// forward mode routine called by CppAD
|
||||
virtual bool forward(
|
||||
// lowest order Taylor coefficient we are evaluating
|
||||
size_t p ,
|
||||
// highest order Taylor coefficient we are evaluating
|
||||
size_t q ,
|
||||
// which components of x are variables
|
||||
const CppAD::vector<bool>& vx ,
|
||||
// which components of y are variables
|
||||
CppAD::vector<bool>& vy ,
|
||||
// tx [ j * (q+1) + k ] is x_j^k
|
||||
const CppAD::vector<scalar>& tx ,
|
||||
// ty [ i * (q+1) + k ] is y_i^k
|
||||
CppAD::vector<scalar>& ty
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
size_t ny = nr * nr;
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 1 + ny;
|
||||
# endif
|
||||
assert( vx.size() == 0 || nx == vx.size() );
|
||||
assert( vx.size() == 0 || ny == vy.size() );
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
//
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_arg_ and f_result_ are large enough
|
||||
assert( f_arg_.size() == f_result_.size() );
|
||||
if( f_arg_.size() < n_order )
|
||||
{ f_arg_.resize(n_order);
|
||||
f_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ f_arg_[k].resize(nr, nr);
|
||||
f_result_[k].resize(nr, nr);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_arg_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // unpack arg values for this order
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// result for each order
|
||||
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
|
||||
//
|
||||
f_result_[0] = f_arg_[0].inverse();
|
||||
for(size_t k = 1; k < n_order; k++)
|
||||
{ // initialize sum
|
||||
matrix f_sum = matrix::Zero(nr, nr);
|
||||
// compute sum
|
||||
for(size_t ell = 1; ell <= k; ell++)
|
||||
f_sum -= f_arg_[ell] * f_result_[k-ell];
|
||||
// result_[k] = arg_[0]^{-1} * sum_
|
||||
f_result_[k] = f_result_[0] * f_sum;
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// pack result_ into ty
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ for(size_t i = 0; i < ny; i++)
|
||||
ty[ i * n_order + k ] = f_result_[k].data()[i];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// check if we are computing vy
|
||||
if( vx.size() == 0 )
|
||||
return true;
|
||||
// ------------------------------------------------------------------
|
||||
// This is a very dumb algorithm that over estimates which
|
||||
// elements of the inverse are variables (which is not efficient).
|
||||
bool var = false;
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
var |= vx[1 + i];
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
vy[i] = var;
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead reverse$$
|
||||
$srccode%cpp% */
|
||||
// reverse mode routine called by CppAD
|
||||
virtual bool reverse(
|
||||
// highest order Taylor coefficient that we are computing derivative of
|
||||
size_t q ,
|
||||
// forward mode Taylor coefficients for x variables
|
||||
const CppAD::vector<double>& tx ,
|
||||
// forward mode Taylor coefficients for y variables
|
||||
const CppAD::vector<double>& ty ,
|
||||
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
|
||||
CppAD::vector<double>& px ,
|
||||
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
|
||||
const CppAD::vector<double>& py
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
size_t ny = nr * nr;
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 1 + ny;
|
||||
# endif
|
||||
//
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
assert( px.size() == tx.size() );
|
||||
assert( py.size() == ty.size() );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_arg_ is large enough
|
||||
assert( f_arg_.size() == f_result_.size() );
|
||||
// must have previous run forward with order >= n_order
|
||||
assert( f_arg_.size() >= n_order );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure r_arg_, r_result_ are large enough
|
||||
assert( r_arg_.size() == r_result_.size() );
|
||||
if( r_arg_.size() < n_order )
|
||||
{ r_arg_.resize(n_order);
|
||||
r_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ r_arg_[k].resize(nr, nr);
|
||||
r_result_[k].resize(nr, nr);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_arg_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // unpack arg values for this order
|
||||
for(size_t i = 0; i < ny; i++)
|
||||
f_arg_[k].data()[i] = tx[ (1 + i) * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack py into r_result_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ for(size_t i = 0; i < ny; i++)
|
||||
r_result_[k].data()[i] = py[ i * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// initialize r_arg_ as zero
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
r_arg_[k] = matrix::Zero(nr, nr);
|
||||
// -------------------------------------------------------------------
|
||||
// matrix reverse mode calculation
|
||||
//
|
||||
for(size_t k1 = n_order; k1 > 1; k1--)
|
||||
{ size_t k = k1 - 1;
|
||||
// bar{R}_0 = bar{R}_0 + bar{R}_k (A_0 R_k)^T
|
||||
r_result_[0] +=
|
||||
r_result_[k] * f_result_[k].transpose() * f_arg_[0].transpose();
|
||||
//
|
||||
for(size_t ell = 1; ell <= k; ell++)
|
||||
{ // bar{A}_l = bar{A}_l - R_0^T bar{R}_k R_{k-l}^T
|
||||
r_arg_[ell] -= f_result_[0].transpose()
|
||||
* r_result_[k] * f_result_[k-ell].transpose();
|
||||
// bar{R}_{k-l} = bar{R}_{k-1} - (R_0 A_l)^T bar{R}_k
|
||||
r_result_[k-ell] -= f_arg_[ell].transpose()
|
||||
* f_result_[0].transpose() * r_result_[k];
|
||||
}
|
||||
}
|
||||
r_arg_[0] -=
|
||||
f_result_[0].transpose() * r_result_[0] * f_result_[0].transpose();
|
||||
// -------------------------------------------------------------------
|
||||
// pack r_arg into px
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ for(size_t i = 0; i < ny; i++)
|
||||
px[ (1 + i) * n_order + k ] = r_arg_[k].data()[i];
|
||||
}
|
||||
//
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head End Class Definition$$
|
||||
$srccode%cpp% */
|
||||
}; // End of atomic_eigen_mat_inv class
|
||||
|
||||
} // END_EMPTY_NAMESPACE
|
||||
/* %$$
|
||||
$$ $comment end nospell$$
|
||||
$end
|
||||
*/
|
||||
|
||||
|
||||
# endif
|
||||
@@ -0,0 +1,656 @@
|
||||
# ifndef CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
|
||||
# define CPPAD_EXAMPLE_EIGEN_MAT_MUL_HPP
|
||||
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
|
||||
/*
|
||||
$begin atomic_eigen_mat_mul.hpp$$
|
||||
$spell
|
||||
Eigen
|
||||
Taylor
|
||||
nr
|
||||
nc
|
||||
$$
|
||||
|
||||
$section Atomic Eigen Matrix Multiply Class$$
|
||||
|
||||
$head See Also$$
|
||||
$cref atomic_mat_mul.hpp$$
|
||||
|
||||
$head Purpose$$
|
||||
Construct an atomic operation that computes the matrix product,
|
||||
$latex R = A \times B$$
|
||||
for any positive integers $latex r$$, $latex m$$, $latex c$$,
|
||||
and any $latex A \in \B{R}^{r \times m}$$,
|
||||
$latex B \in \B{R}^{m \times c}$$.
|
||||
|
||||
$head Matrix Dimensions$$
|
||||
This example puts the matrix dimensions in the atomic function arguments,
|
||||
instead of the $cref/constructor/atomic_ctor/$$, so that they can
|
||||
be different for different calls to the atomic function.
|
||||
These dimensions are:
|
||||
$table
|
||||
$icode nr_left$$
|
||||
$cnext number of rows in the left matrix; i.e, $latex r$$ $rend
|
||||
$icode n_middle$$
|
||||
$cnext rows in the left matrix and columns in right; i.e, $latex m$$ $rend
|
||||
$icode nc_right$$
|
||||
$cnext number of columns in the right matrix; i.e., $latex c$$
|
||||
$tend
|
||||
|
||||
$head Theory$$
|
||||
|
||||
$subhead Forward$$
|
||||
For $latex k = 0 , \ldots $$, the $th k$$ order Taylor coefficient
|
||||
$latex R_k$$ is given by
|
||||
$latex \[
|
||||
R_k = \sum_{\ell = 0}^{k} A_\ell B_{k-\ell}
|
||||
\] $$
|
||||
|
||||
$subhead Product of Two Matrices$$
|
||||
Suppose $latex \bar{E}$$ is the derivative of the
|
||||
scalar value function $latex s(E)$$ with respect to $latex E$$; i.e.,
|
||||
$latex \[
|
||||
\bar{E}_{i,j} = \frac{ \partial s } { \partial E_{i,j} }
|
||||
\] $$
|
||||
Also suppose that $latex t$$ is a scalar valued argument and
|
||||
$latex \[
|
||||
E(t) = C(t) D(t)
|
||||
\] $$
|
||||
It follows that
|
||||
$latex \[
|
||||
E'(t) = C'(t) D(t) + C(t) D'(t)
|
||||
\] $$
|
||||
|
||||
$latex \[
|
||||
(s \circ E)'(t)
|
||||
=
|
||||
\R{tr} [ \bar{E}^\R{T} E'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
=
|
||||
\R{tr} [ \bar{E}^\R{T} C'(t) D(t) ] +
|
||||
\R{tr} [ \bar{E}^\R{T} C(t) D'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
=
|
||||
\R{tr} [ D(t) \bar{E}^\R{T} C'(t) ] +
|
||||
\R{tr} [ \bar{E}^\R{T} C(t) D'(t) ]
|
||||
\] $$
|
||||
$latex \[
|
||||
\bar{C} = \bar{E} D^\R{T} \W{,}
|
||||
\bar{D} = C^\R{T} \bar{E}
|
||||
\] $$
|
||||
|
||||
$subhead Reverse$$
|
||||
Reverse mode eliminates $latex R_k$$ as follows:
|
||||
for $latex \ell = 0, \ldots , k-1$$,
|
||||
$latex \[
|
||||
\bar{A}_\ell = \bar{A}_\ell + \bar{R}_k B_{k-\ell}^\R{T}
|
||||
\] $$
|
||||
$latex \[
|
||||
\bar{B}_{k-\ell} = \bar{B}_{k-\ell} + A_\ell^\R{T} \bar{R}_k
|
||||
\] $$
|
||||
|
||||
|
||||
$nospell
|
||||
|
||||
$head Start Class Definition$$
|
||||
$srccode%cpp% */
|
||||
# include <cppad/cppad.hpp>
|
||||
# include <Eigen/Core>
|
||||
|
||||
|
||||
/* %$$
|
||||
$head Public$$
|
||||
|
||||
$subhead Types$$
|
||||
$srccode%cpp% */
|
||||
namespace { // BEGIN_EMPTY_NAMESPACE
|
||||
|
||||
template <class Base>
|
||||
class atomic_eigen_mat_mul : public CppAD::atomic_base<Base> {
|
||||
public:
|
||||
// -----------------------------------------------------------
|
||||
// type of elements during calculation of derivatives
|
||||
typedef Base scalar;
|
||||
// type of elements during taping
|
||||
typedef CppAD::AD<scalar> ad_scalar;
|
||||
// type of matrix during calculation of derivatives
|
||||
typedef Eigen::Matrix<
|
||||
scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> matrix;
|
||||
// type of matrix during taping
|
||||
typedef Eigen::Matrix<
|
||||
ad_scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ad_matrix;
|
||||
/* %$$
|
||||
$subhead Constructor$$
|
||||
$srccode%cpp% */
|
||||
// constructor
|
||||
atomic_eigen_mat_mul(void) : CppAD::atomic_base<Base>(
|
||||
"atom_eigen_mat_mul" ,
|
||||
CppAD::atomic_base<Base>::set_sparsity_enum
|
||||
)
|
||||
{ }
|
||||
/* %$$
|
||||
$subhead op$$
|
||||
$srccode%cpp% */
|
||||
// use atomic operation to multiply two AD matrices
|
||||
ad_matrix op(
|
||||
const ad_matrix& left ,
|
||||
const ad_matrix& right )
|
||||
{ size_t nr_left = size_t( left.rows() );
|
||||
size_t n_middle = size_t( left.cols() );
|
||||
size_t nc_right = size_t( right.cols() );
|
||||
assert( n_middle == size_t( right.rows() ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
size_t n_left = nr_left * n_middle;
|
||||
size_t n_right = n_middle * nc_right;
|
||||
size_t n_result = nr_left * nc_right;
|
||||
//
|
||||
assert( 3 + n_left + n_right == nx );
|
||||
assert( n_result == ny );
|
||||
// -----------------------------------------------------------------
|
||||
// packed version of left and right
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_arg(nx);
|
||||
//
|
||||
packed_arg[0] = ad_scalar( nr_left );
|
||||
packed_arg[1] = ad_scalar( n_middle );
|
||||
packed_arg[2] = ad_scalar( nc_right );
|
||||
for(size_t i = 0; i < n_left; i++)
|
||||
packed_arg[3 + i] = left.data()[i];
|
||||
for(size_t i = 0; i < n_right; i++)
|
||||
packed_arg[ 3 + n_left + i ] = right.data()[i];
|
||||
// ------------------------------------------------------------------
|
||||
// Packed version of result = left * right.
|
||||
// This as an atomic_base funciton call that CppAD uses
|
||||
// to store the atomic operation on the tape.
|
||||
CPPAD_TESTVECTOR(ad_scalar) packed_result(ny);
|
||||
(*this)(packed_arg, packed_result);
|
||||
// ------------------------------------------------------------------
|
||||
// unpack result matrix
|
||||
ad_matrix result(nr_left, nc_right);
|
||||
for(size_t i = 0; i < n_result; i++)
|
||||
result.data()[i] = packed_result[ i ];
|
||||
//
|
||||
return result;
|
||||
}
|
||||
/* %$$
|
||||
$head Private$$
|
||||
|
||||
$subhead Variables$$
|
||||
$srccode%cpp% */
|
||||
private:
|
||||
// -------------------------------------------------------------
|
||||
// one forward mode vector of matrices for left, right, and result
|
||||
CppAD::vector<matrix> f_left_, f_right_, f_result_;
|
||||
// one reverse mode vector of matrices for left, right, and result
|
||||
CppAD::vector<matrix> r_left_, r_right_, r_result_;
|
||||
// -------------------------------------------------------------
|
||||
/* %$$
|
||||
$subhead forward$$
|
||||
$srccode%cpp% */
|
||||
// forward mode routine called by CppAD
|
||||
virtual bool forward(
|
||||
// lowest order Taylor coefficient we are evaluating
|
||||
size_t p ,
|
||||
// highest order Taylor coefficient we are evaluating
|
||||
size_t q ,
|
||||
// which components of x are variables
|
||||
const CppAD::vector<bool>& vx ,
|
||||
// which components of y are variables
|
||||
CppAD::vector<bool>& vy ,
|
||||
// tx [ 3 + j * (q+1) + k ] is x_j^k
|
||||
const CppAD::vector<scalar>& tx ,
|
||||
// ty [ i * (q+1) + k ] is y_i^k
|
||||
CppAD::vector<scalar>& ty
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( vx.size() == 0 || nx == vx.size() );
|
||||
assert( vx.size() == 0 || ny == vy.size() );
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
//
|
||||
size_t n_left = nr_left * n_middle;
|
||||
size_t n_right = n_middle * nc_right;
|
||||
size_t n_result = nr_left * nc_right;
|
||||
assert( 3 + n_left + n_right == nx );
|
||||
assert( n_result == ny );
|
||||
//
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_left_, f_right_, and f_result_ are large enough
|
||||
assert( f_left_.size() == f_right_.size() );
|
||||
assert( f_left_.size() == f_result_.size() );
|
||||
if( f_left_.size() < n_order )
|
||||
{ f_left_.resize(n_order);
|
||||
f_right_.resize(n_order);
|
||||
f_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ f_left_[k].resize(nr_left, n_middle);
|
||||
f_right_[k].resize(n_middle, nc_right);
|
||||
f_result_[k].resize(nr_left, nc_right);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_left and f_right
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // unpack left values for this order
|
||||
for(size_t i = 0; i < n_left; i++)
|
||||
f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
|
||||
//
|
||||
// unpack right values for this order
|
||||
for(size_t i = 0; i < n_right; i++)
|
||||
f_right_[k].data()[i] = tx[ ( 3 + n_left + i) * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// result for each order
|
||||
// (we could avoid recalculting f_result_[k] for k=0,...,p-1)
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // result[k] = sum_ell left[ell] * right[k-ell]
|
||||
f_result_[k] = matrix::Zero(nr_left, nc_right);
|
||||
for(size_t ell = 0; ell <= k; ell++)
|
||||
f_result_[k] += f_left_[ell] * f_right_[k-ell];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// pack result_ into ty
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ for(size_t i = 0; i < n_result; i++)
|
||||
ty[ i * n_order + k ] = f_result_[k].data()[i];
|
||||
}
|
||||
// ------------------------------------------------------------------
|
||||
// check if we are computing vy
|
||||
if( vx.size() == 0 )
|
||||
return true;
|
||||
// ------------------------------------------------------------------
|
||||
// compute variable information for y; i.e., vy
|
||||
// (note that the constant zero times a variable is a constant)
|
||||
scalar zero(0.0);
|
||||
assert( n_order == 1 );
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ bool var = false;
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ // left information
|
||||
size_t index = 3 + i * n_middle + ell;
|
||||
bool var_left = vx[index];
|
||||
bool nz_left = var_left | (f_left_[0](i, ell) != zero);
|
||||
// right information
|
||||
index = 3 + n_left + ell * nc_right + j;
|
||||
bool var_right = vx[index];
|
||||
bool nz_right = var_right | (f_right_[0](ell, j) != zero);
|
||||
// effect of result
|
||||
var |= var_left & nz_right;
|
||||
var |= nz_left & var_right;
|
||||
}
|
||||
size_t index = i * nc_right + j;
|
||||
vy[index] = var;
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead reverse$$
|
||||
$srccode%cpp% */
|
||||
// reverse mode routine called by CppAD
|
||||
virtual bool reverse(
|
||||
// highest order Taylor coefficient that we are computing derivative of
|
||||
size_t q ,
|
||||
// forward mode Taylor coefficients for x variables
|
||||
const CppAD::vector<double>& tx ,
|
||||
// forward mode Taylor coefficients for y variables
|
||||
const CppAD::vector<double>& ty ,
|
||||
// upon return, derivative of G[ F[ {x_j^k} ] ] w.r.t {x_j^k}
|
||||
CppAD::vector<double>& px ,
|
||||
// derivative of G[ {y_i^k} ] w.r.t. {y_i^k}
|
||||
const CppAD::vector<double>& py
|
||||
)
|
||||
{ size_t n_order = q + 1;
|
||||
size_t nr_left = size_t( CppAD::Integer( tx[ 0 * n_order + 0 ] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( tx[ 1 * n_order + 0 ] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( tx[ 2 * n_order + 0 ] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
assert( px.size() == tx.size() );
|
||||
assert( py.size() == ty.size() );
|
||||
//
|
||||
size_t n_left = nr_left * n_middle;
|
||||
size_t n_right = n_middle * nc_right;
|
||||
size_t n_result = nr_left * nc_right;
|
||||
assert( 3 + n_left + n_right == nx );
|
||||
assert( n_result == ny );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure f_left_, f_right_ are large enough
|
||||
assert( f_left_.size() == f_right_.size() );
|
||||
assert( f_left_.size() == f_result_.size() );
|
||||
// must have previous run forward with order >= n_order
|
||||
assert( f_left_.size() >= n_order );
|
||||
// -------------------------------------------------------------------
|
||||
// make sure r_left_, r_right_, and r_result_ are large enough
|
||||
assert( r_left_.size() == r_right_.size() );
|
||||
assert( r_left_.size() == r_result_.size() );
|
||||
if( r_left_.size() < n_order )
|
||||
{ r_left_.resize(n_order);
|
||||
r_right_.resize(n_order);
|
||||
r_result_.resize(n_order);
|
||||
//
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ r_left_[k].resize(nr_left, n_middle);
|
||||
r_right_[k].resize(n_middle, nc_right);
|
||||
r_result_[k].resize(nr_left, nc_right);
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack tx into f_left and f_right
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // unpack left values for this order
|
||||
for(size_t i = 0; i < n_left; i++)
|
||||
f_left_[k].data()[i] = tx[ (3 + i) * n_order + k ];
|
||||
//
|
||||
// unpack right values for this order
|
||||
for(size_t i = 0; i < n_right; i++)
|
||||
f_right_[k].data()[i] = tx[ (3 + n_left + i) * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// unpack py into r_result_
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ for(size_t i = 0; i < n_result; i++)
|
||||
r_result_[k].data()[i] = py[ i * n_order + k ];
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// initialize r_left_ and r_right_ as zero
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ r_left_[k] = matrix::Zero(nr_left, n_middle);
|
||||
r_right_[k] = matrix::Zero(n_middle, nc_right);
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// matrix reverse mode calculation
|
||||
for(size_t k1 = n_order; k1 > 0; k1--)
|
||||
{ size_t k = k1 - 1;
|
||||
for(size_t ell = 0; ell <= k; ell++)
|
||||
{ // nr x nm = nr x nc * nc * nm
|
||||
r_left_[ell] += r_result_[k] * f_right_[k-ell].transpose();
|
||||
// nm x nc = nm x nr * nr * nc
|
||||
r_right_[k-ell] += f_left_[ell].transpose() * r_result_[k];
|
||||
}
|
||||
}
|
||||
// -------------------------------------------------------------------
|
||||
// pack r_left and r_right int px
|
||||
for(size_t k = 0; k < n_order; k++)
|
||||
{ // dimensions are integer constants
|
||||
px[ 0 * n_order + k ] = 0.0;
|
||||
px[ 1 * n_order + k ] = 0.0;
|
||||
px[ 2 * n_order + k ] = 0.0;
|
||||
//
|
||||
// pack left values for this order
|
||||
for(size_t i = 0; i < n_left; i++)
|
||||
px[ (3 + i) * n_order + k ] = r_left_[k].data()[i];
|
||||
//
|
||||
// pack right values for this order
|
||||
for(size_t i = 0; i < n_right; i++)
|
||||
px[ (3 + i + n_left) * n_order + k] = r_right_[k].data()[i];
|
||||
}
|
||||
//
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead for_sparse_jac$$
|
||||
$srccode%cpp% */
|
||||
// forward Jacobian sparsity routine called by CppAD
|
||||
virtual bool for_sparse_jac(
|
||||
// number of columns in the matrix R
|
||||
size_t q ,
|
||||
// sparsity pattern for the matrix R
|
||||
const CppAD::vector< std::set<size_t> >& r ,
|
||||
// sparsity pattern for the matrix S = f'(x) * R
|
||||
CppAD::vector< std::set<size_t> >& s ,
|
||||
const CppAD::vector<Base>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( nx == r.size() );
|
||||
assert( ny == s.size() );
|
||||
//
|
||||
size_t n_left = nr_left * n_middle;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ // pack index for entry (i, j) in result
|
||||
size_t i_result = i * nc_right + j;
|
||||
s[i_result].clear();
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ // pack index for entry (i, ell) in left
|
||||
size_t i_left = 3 + i * n_middle + ell;
|
||||
// pack index for entry (ell, j) in right
|
||||
size_t i_right = 3 + n_left + ell * nc_right + j;
|
||||
// check if result of for this product is alwasy zero
|
||||
// note that x is nan for commponents that are variables
|
||||
bool zero = x[i_left] == Base(0.0) || x[i_right] == Base(0);
|
||||
if( ! zero )
|
||||
{ s[i_result] =
|
||||
CppAD::set_union(s[i_result], r[i_left] );
|
||||
s[i_result] =
|
||||
CppAD::set_union(s[i_result], r[i_right] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead rev_sparse_jac$$
|
||||
$srccode%cpp% */
|
||||
// reverse Jacobian sparsity routine called by CppAD
|
||||
virtual bool rev_sparse_jac(
|
||||
// number of columns in the matrix R^T
|
||||
size_t q ,
|
||||
// sparsity pattern for the matrix R^T
|
||||
const CppAD::vector< std::set<size_t> >& rt ,
|
||||
// sparsoity pattern for the matrix S^T = f'(x)^T * R^T
|
||||
CppAD::vector< std::set<size_t> >& st ,
|
||||
const CppAD::vector<Base>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( nx == st.size() );
|
||||
assert( ny == rt.size() );
|
||||
//
|
||||
// initialize S^T as empty
|
||||
for(size_t i = 0; i < nx; i++)
|
||||
st[i].clear();
|
||||
|
||||
// sparsity for S(x)^T = f'(x)^T * R^T
|
||||
size_t n_left = nr_left * n_middle;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ // pack index for entry (i, j) in result
|
||||
size_t i_result = i * nc_right + j;
|
||||
st[i_result].clear();
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ // pack index for entry (i, ell) in left
|
||||
size_t i_left = 3 + i * n_middle + ell;
|
||||
// pack index for entry (ell, j) in right
|
||||
size_t i_right = 3 + n_left + ell * nc_right + j;
|
||||
//
|
||||
st[i_left] = CppAD::set_union(st[i_left], rt[i_result]);
|
||||
st[i_right] = CppAD::set_union(st[i_right], rt[i_result]);
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead for_sparse_hes$$
|
||||
$srccode%cpp% */
|
||||
virtual bool for_sparse_hes(
|
||||
// which components of x are variables for this call
|
||||
const CppAD::vector<bool>& vx,
|
||||
// sparsity pattern for the diagonal of R
|
||||
const CppAD::vector<bool>& r ,
|
||||
// sparsity pattern for the vector S
|
||||
const CppAD::vector<bool>& s ,
|
||||
// sparsity patternfor the Hessian H(x)
|
||||
CppAD::vector< std::set<size_t> >& h ,
|
||||
const CppAD::vector<Base>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( vx.size() == nx );
|
||||
assert( r.size() == nx );
|
||||
assert( s.size() == ny );
|
||||
assert( h.size() == nx );
|
||||
//
|
||||
// initilize h as empty
|
||||
for(size_t i = 0; i < nx; i++)
|
||||
h[i].clear();
|
||||
//
|
||||
size_t n_left = nr_left * n_middle;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ // pack index for entry (i, j) in result
|
||||
size_t i_result = i * nc_right + j;
|
||||
if( s[i_result] )
|
||||
{ for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ // pack index for entry (i, ell) in left
|
||||
size_t i_left = 3 + i * n_middle + ell;
|
||||
// pack index for entry (ell, j) in right
|
||||
size_t i_right = 3 + n_left + ell * nc_right + j;
|
||||
if( r[i_left] & r[i_right] )
|
||||
{ h[i_left].insert(i_right);
|
||||
h[i_right].insert(i_left);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$subhead rev_sparse_hes$$
|
||||
$srccode%cpp% */
|
||||
// reverse Hessian sparsity routine called by CppAD
|
||||
virtual bool rev_sparse_hes(
|
||||
// which components of x are variables for this call
|
||||
const CppAD::vector<bool>& vx,
|
||||
// sparsity pattern for S(x) = g'[f(x)]
|
||||
const CppAD::vector<bool>& s ,
|
||||
// sparsity pattern for d/dx g[f(x)] = S(x) * f'(x)
|
||||
CppAD::vector<bool>& t ,
|
||||
// number of columns in R, U(x), and V(x)
|
||||
size_t q ,
|
||||
// sparsity pattern for R
|
||||
const CppAD::vector< std::set<size_t> >& r ,
|
||||
// sparsity pattern for U(x) = g^{(2)} [ f(x) ] * f'(x) * R
|
||||
const CppAD::vector< std::set<size_t> >& u ,
|
||||
// sparsity pattern for
|
||||
// V(x) = f'(x)^T * U(x) + sum_{i=0}^{m-1} S_i(x) f_i^{(2)} (x) * R
|
||||
CppAD::vector< std::set<size_t> >& v ,
|
||||
// parameters as integers
|
||||
const CppAD::vector<Base>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
//
|
||||
assert( vx.size() == nx );
|
||||
assert( s.size() == ny );
|
||||
assert( t.size() == nx );
|
||||
assert( r.size() == nx );
|
||||
assert( v.size() == nx );
|
||||
//
|
||||
// initilaize return sparsity patterns as false
|
||||
for(size_t j = 0; j < nx; j++)
|
||||
{ t[j] = false;
|
||||
v[j].clear();
|
||||
}
|
||||
//
|
||||
size_t n_left = nr_left * n_middle;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ // pack index for entry (i, j) in result
|
||||
size_t i_result = i * nc_right + j;
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ // pack index for entry (i, ell) in left
|
||||
size_t i_left = 3 + i * n_middle + ell;
|
||||
// pack index for entry (ell, j) in right
|
||||
size_t i_right = 3 + n_left + ell * nc_right + j;
|
||||
//
|
||||
// back propagate T(x) = S(x) * f'(x).
|
||||
t[i_left] |= bool( s[i_result] );
|
||||
t[i_right] |= bool( s[i_result] );
|
||||
//
|
||||
// V(x) = f'(x)^T * U(x) + sum_i S_i(x) * f_i''(x) * R
|
||||
// U(x) = g''[ f(x) ] * f'(x) * R
|
||||
// S_i(x) = g_i'[ f(x) ]
|
||||
//
|
||||
// back propagate f'(x)^T * U(x)
|
||||
v[i_left] = CppAD::set_union(v[i_left], u[i_result] );
|
||||
v[i_right] = CppAD::set_union(v[i_right], u[i_result] );
|
||||
//
|
||||
// back propagate S_i(x) * f_i''(x) * R
|
||||
// (here is where we use vx to check for cross terms)
|
||||
if( s[i_result] & vx[i_left] & vx[i_right] )
|
||||
{ v[i_left] = CppAD::set_union(v[i_left], r[i_right] );
|
||||
v[i_right] = CppAD::set_union(v[i_right], r[i_left] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head End Class Definition$$
|
||||
$srccode%cpp% */
|
||||
}; // End of atomic_eigen_mat_mul class
|
||||
|
||||
} // END_EMPTY_NAMESPACE
|
||||
/* %$$
|
||||
$$ $comment end nospell$$
|
||||
$end
|
||||
*/
|
||||
|
||||
|
||||
# endif
|
||||
+678
@@ -0,0 +1,678 @@
|
||||
// $Id$
|
||||
# ifndef CPPAD_EXAMPLE_MAT_MUL_HPP
|
||||
# define CPPAD_EXAMPLE_MAT_MUL_HPP
|
||||
|
||||
/* --------------------------------------------------------------------------
|
||||
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
|
||||
|
||||
CppAD is distributed under multiple licenses. This distribution is under
|
||||
the terms of the
|
||||
Eclipse Public License Version 1.0.
|
||||
|
||||
A copy of this license is included in the COPYING file of this distribution.
|
||||
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
||||
-------------------------------------------------------------------------- */
|
||||
|
||||
/*
|
||||
$begin atomic_mat_mul.hpp$$
|
||||
$spell
|
||||
Taylor
|
||||
ty
|
||||
px
|
||||
CppAD
|
||||
jac
|
||||
hes
|
||||
nr
|
||||
nc
|
||||
$$
|
||||
|
||||
$section Matrix Multiply as an Atomic Operation$$
|
||||
|
||||
$head See Also$$
|
||||
$cref atomic_eigen_mat_mul.hpp$$
|
||||
|
||||
$head Matrix Dimensions$$
|
||||
This example puts the matrix dimensions in the atomic function arguments,
|
||||
instead of the $cref/constructor/atomic_ctor/$$, so that they can
|
||||
be different for different calls to the atomic function.
|
||||
These dimensions are:
|
||||
$table
|
||||
$icode nr_left$$ $cnext number of rows in the left matrix $rend
|
||||
$icode n_middle$$ $cnext rows in the left matrix and columns in right $rend
|
||||
$icode nc_right$$ $cnext number of columns in the right matrix
|
||||
$tend
|
||||
|
||||
$head Start Class Definition$$
|
||||
$srccode%cpp% */
|
||||
# include <cppad/cppad.hpp>
|
||||
namespace { // Begin empty namespace
|
||||
using CppAD::vector;
|
||||
//
|
||||
using CppAD::set_union;
|
||||
//
|
||||
// matrix result = left * right
|
||||
class atomic_mat_mul : public CppAD::atomic_base<double> {
|
||||
/* %$$
|
||||
$head Constructor$$
|
||||
$srccode%cpp% */
|
||||
public:
|
||||
// ---------------------------------------------------------------------
|
||||
// constructor
|
||||
atomic_mat_mul(void) : CppAD::atomic_base<double>("mat_mul")
|
||||
{ }
|
||||
private:
|
||||
/* %$$
|
||||
$head Left Operand Element Index$$
|
||||
Index in the Taylor coefficient matrix $icode tx$$ of a left matrix element.
|
||||
$srccode%cpp% */
|
||||
size_t left(
|
||||
size_t i , // left matrix row index
|
||||
size_t j , // left matrix column index
|
||||
size_t k , // Taylor coeffocient order
|
||||
size_t nk , // number of Taylor coefficients in tx
|
||||
size_t nr_left , // rows in left matrix
|
||||
size_t n_middle , // rows in left and columns in right
|
||||
size_t nc_right ) // columns in right matrix
|
||||
{ assert( i < nr_left );
|
||||
assert( j < n_middle );
|
||||
return (3 + i * n_middle + j) * nk + k;
|
||||
}
|
||||
/* %$$
|
||||
$head Right Operand Element Index$$
|
||||
Index in the Taylor coefficient matrix $icode tx$$ of a right matrix element.
|
||||
$srccode%cpp% */
|
||||
size_t right(
|
||||
size_t i , // right matrix row index
|
||||
size_t j , // right matrix column index
|
||||
size_t k , // Taylor coeffocient order
|
||||
size_t nk , // number of Taylor coefficients in tx
|
||||
size_t nr_left , // rows in left matrix
|
||||
size_t n_middle , // rows in left and columns in right
|
||||
size_t nc_right ) // columns in right matrix
|
||||
{ assert( i < n_middle );
|
||||
assert( j < nc_right );
|
||||
size_t offset = 3 + nr_left * n_middle;
|
||||
return (offset + i * nc_right + j) * nk + k;
|
||||
}
|
||||
/* %$$
|
||||
$head Result Element Index$$
|
||||
Index in the Taylor coefficient matrix $icode ty$$ of a result matrix element.
|
||||
$srccode%cpp% */
|
||||
size_t result(
|
||||
size_t i , // result matrix row index
|
||||
size_t j , // result matrix column index
|
||||
size_t k , // Taylor coeffocient order
|
||||
size_t nk , // number of Taylor coefficients in ty
|
||||
size_t nr_left , // rows in left matrix
|
||||
size_t n_middle , // rows in left and columns in right
|
||||
size_t nc_right ) // columns in right matrix
|
||||
{ assert( i < nr_left );
|
||||
assert( j < nc_right );
|
||||
return (i * nc_right + j) * nk + k;
|
||||
}
|
||||
/* %$$
|
||||
$head Forward Matrix Multiply$$
|
||||
Forward mode multiply Taylor coefficients in $icode tx$$ and sum into
|
||||
$icode ty$$ (for one pair of left and right orders)
|
||||
$srccode%cpp% */
|
||||
void forward_multiply(
|
||||
size_t k_left , // order for left coefficients
|
||||
size_t k_right , // order for right coefficients
|
||||
const vector<double>& tx , // domain space Taylor coefficients
|
||||
vector<double>& ty , // range space Taylor coefficients
|
||||
size_t nr_left , // rows in left matrix
|
||||
size_t n_middle , // rows in left and columns in right
|
||||
size_t nc_right ) // columns in right matrix
|
||||
{
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t nk = tx.size() / nx;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
assert( nk == ty.size() / ny );
|
||||
# endif
|
||||
//
|
||||
size_t k_result = k_left + k_right;
|
||||
assert( k_result < nk );
|
||||
//
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ double sum = 0.0;
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k_left, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k_right, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
sum += tx[i_left] * tx[i_right];
|
||||
}
|
||||
size_t i_result = result(
|
||||
i, j, k_result, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
ty[i_result] += sum;
|
||||
}
|
||||
}
|
||||
}
|
||||
/* %$$
|
||||
$head Reverse Matrix Multiply$$
|
||||
Reverse mode partials of Taylor coefficients and sum into $icode px$$
|
||||
(for one pair of left and right orders)
|
||||
$srccode%cpp% */
|
||||
void reverse_multiply(
|
||||
size_t k_left , // order for left coefficients
|
||||
size_t k_right , // order for right coefficients
|
||||
const vector<double>& tx , // domain space Taylor coefficients
|
||||
const vector<double>& ty , // range space Taylor coefficients
|
||||
vector<double>& px , // partials w.r.t. tx
|
||||
const vector<double>& py , // partials w.r.t. ty
|
||||
size_t nr_left , // rows in left matrix
|
||||
size_t n_middle , // rows in left and columns in right
|
||||
size_t nc_right ) // columns in right matrix
|
||||
{
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t nk = tx.size() / nx;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
assert( nk == ty.size() / ny );
|
||||
# endif
|
||||
assert( tx.size() == px.size() );
|
||||
assert( ty.size() == py.size() );
|
||||
//
|
||||
size_t k_result = k_left + k_right;
|
||||
assert( k_result < nk );
|
||||
//
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k_result, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k_left, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k_right, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
// sum += tx[i_left] * tx[i_right];
|
||||
px[i_left] += tx[i_right] * py[i_result];
|
||||
px[i_right] += tx[i_left] * py[i_result];
|
||||
}
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
/* %$$
|
||||
$head forward$$
|
||||
Routine called by CppAD during $cref Forward$$ mode.
|
||||
$srccode%cpp% */
|
||||
virtual bool forward(
|
||||
size_t q ,
|
||||
size_t p ,
|
||||
const vector<bool>& vx ,
|
||||
vector<bool>& vy ,
|
||||
const vector<double>& tx ,
|
||||
vector<double>& ty
|
||||
)
|
||||
{ size_t n_order = p + 1;
|
||||
size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
|
||||
size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
|
||||
size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( vx.size() == 0 || nx == vx.size() );
|
||||
assert( vx.size() == 0 || ny == vy.size() );
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
size_t i, j, ell;
|
||||
|
||||
// check if we are computing vy information
|
||||
if( vx.size() > 0 )
|
||||
{ size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ bool var = false;
|
||||
for(ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
bool nz_left = vx[i_left] |(tx[i_left] != 0.);
|
||||
bool nz_right = vx[i_right]|(tx[i_right] != 0.);
|
||||
// if not multiplying by the constant zero
|
||||
if( nz_left & nz_right )
|
||||
var |= bool(vx[i_left]) | bool(vx[i_right]);
|
||||
}
|
||||
size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
vy[i_result] = var;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// initialize result as zero
|
||||
size_t k;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ for(k = q; k <= p; k++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, n_order, nr_left, n_middle, nc_right
|
||||
);
|
||||
ty[i_result] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
for(k = q; k <= p; k++)
|
||||
{ // sum the produces that result in order k
|
||||
for(ell = 0; ell <= k; ell++)
|
||||
forward_multiply(
|
||||
ell, k - ell, tx, ty, nr_left, n_middle, nc_right
|
||||
);
|
||||
}
|
||||
|
||||
// all orders are implented, so always return true
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head reverse$$
|
||||
Routine called by CppAD during $cref Reverse$$ mode.
|
||||
$srccode%cpp% */
|
||||
virtual bool reverse(
|
||||
size_t p ,
|
||||
const vector<double>& tx ,
|
||||
const vector<double>& ty ,
|
||||
vector<double>& px ,
|
||||
const vector<double>& py
|
||||
)
|
||||
{ size_t n_order = p + 1;
|
||||
size_t nr_left = size_t( tx[ 0 * n_order + 0 ] );
|
||||
size_t n_middle = size_t( tx[ 1 * n_order + 0 ] );
|
||||
size_t nc_right = size_t( tx[ 2 * n_order + 0 ] );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( nx * n_order == tx.size() );
|
||||
assert( ny * n_order == ty.size() );
|
||||
assert( px.size() == tx.size() );
|
||||
assert( py.size() == ty.size() );
|
||||
|
||||
// initialize summation
|
||||
for(size_t i = 0; i < px.size(); i++)
|
||||
px[i] = 0.0;
|
||||
|
||||
// number of orders to differentiate
|
||||
size_t k = n_order;
|
||||
while(k--)
|
||||
{ // differentiate the produces that result in order k
|
||||
for(size_t ell = 0; ell <= k; ell++)
|
||||
reverse_multiply(
|
||||
ell, k - ell, tx, ty, px, py, nr_left, n_middle, nc_right
|
||||
);
|
||||
}
|
||||
|
||||
// all orders are implented, so always return true
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head for_sparse_jac$$
|
||||
Routines called by CppAD during $cref ForSparseJac$$.
|
||||
$srccode%cpp% */
|
||||
// boolean sparsity patterns
|
||||
virtual bool for_sparse_jac(
|
||||
size_t q ,
|
||||
const vector<bool>& r ,
|
||||
vector<bool>& s ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( nx == x.size() );
|
||||
assert( nx * q == r.size() );
|
||||
assert( ny * q == s.size() );
|
||||
size_t p;
|
||||
|
||||
// sparsity for S(x) = f'(x) * R
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(p = 0; p < q; p++)
|
||||
s[i_result * q + p] = false;
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(p = 0; p < q; p++)
|
||||
{ // cast avoids Microsoft warning (should not be needed)
|
||||
s[i_result * q + p] |= bool( r[i_left * q + p ] );
|
||||
s[i_result * q + p] |= bool( r[i_right * q + p ] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
// set sparsity patterns
|
||||
virtual bool for_sparse_jac(
|
||||
size_t q ,
|
||||
const vector< std::set<size_t> >& r ,
|
||||
vector< std::set<size_t> >& s ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
# ifndef NDEBUG
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( nx == x.size() );
|
||||
assert( nx == r.size() );
|
||||
assert( ny == s.size() );
|
||||
|
||||
// sparsity for S(x) = f'(x) * R
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(size_t i = 0; i < nr_left; i++)
|
||||
{ for(size_t j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
s[i_result].clear();
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
//
|
||||
s[i_result] = set_union(s[i_result], r[i_left] );
|
||||
s[i_result] = set_union(s[i_result], r[i_right] );
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head rev_sparse_jac$$
|
||||
Routines called by CppAD during $cref RevSparseJac$$.
|
||||
$srccode%cpp% */
|
||||
// boolean sparsity patterns
|
||||
virtual bool rev_sparse_jac(
|
||||
size_t q ,
|
||||
const vector<bool>& rt ,
|
||||
vector<bool>& st ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( nx == x.size() );
|
||||
assert( nx * q == st.size() );
|
||||
assert( ny * q == rt.size() );
|
||||
size_t i, j, p;
|
||||
|
||||
// initialize
|
||||
for(i = 0; i < nx; i++)
|
||||
{ for(p = 0; p < q; p++)
|
||||
st[ i * q + p ] = false;
|
||||
}
|
||||
|
||||
// sparsity for S(x)^T = f'(x)^T * R^T
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(p = 0; p < q; p++)
|
||||
{ st[i_left * q + p] |= bool( rt[i_result * q + p] );
|
||||
st[i_right* q + p] |= bool( rt[i_result * q + p] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
// set sparsity patterns
|
||||
virtual bool rev_sparse_jac(
|
||||
size_t q ,
|
||||
const vector< std::set<size_t> >& rt ,
|
||||
vector< std::set<size_t> >& st ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( nx == x.size() );
|
||||
assert( nx == st.size() );
|
||||
assert( ny == rt.size() );
|
||||
size_t i, j;
|
||||
|
||||
// initialize
|
||||
for(i = 0; i < nx; i++)
|
||||
st[i].clear();
|
||||
|
||||
// sparsity for S(x)^T = f'(x)^T * R^T
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
//
|
||||
st[i_left] = set_union(st[i_left], rt[i_result]);
|
||||
st[i_right] = set_union(st[i_right], rt[i_result]);
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head rev_sparse_hes$$
|
||||
Routines called by $cref RevSparseHes$$.
|
||||
$srccode%cpp% */
|
||||
// set sparsity patterns
|
||||
virtual bool rev_sparse_hes(
|
||||
const vector<bool>& vx,
|
||||
const vector<bool>& s ,
|
||||
vector<bool>& t ,
|
||||
size_t q ,
|
||||
const vector< std::set<size_t> >& r ,
|
||||
const vector< std::set<size_t> >& u ,
|
||||
vector< std::set<size_t> >& v ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( x.size() == nx );
|
||||
assert( vx.size() == nx );
|
||||
assert( t.size() == nx );
|
||||
assert( r.size() == nx );
|
||||
assert( v.size() == nx );
|
||||
assert( s.size() == ny );
|
||||
assert( u.size() == ny );
|
||||
//
|
||||
size_t i, j;
|
||||
//
|
||||
// initilaize sparsity patterns as false
|
||||
for(j = 0; j < nx; j++)
|
||||
{ t[j] = false;
|
||||
v[j].clear();
|
||||
}
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
//
|
||||
// Compute sparsity for T(x) = S(x) * f'(x).
|
||||
// We need not use vx with f'(x) back propagation.
|
||||
t[i_left] |= bool( s[i_result] );
|
||||
t[i_right] |= bool( s[i_result] );
|
||||
|
||||
// V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
|
||||
// U(x) = g''(y) * f'(x) * R
|
||||
// S(x) = g'(y)
|
||||
|
||||
// back propagate f'(x)^T * U(x)
|
||||
// (no need to use vx with f'(x) propogation)
|
||||
v[i_left] = set_union(v[i_left], u[i_result] );
|
||||
v[i_right] = set_union(v[i_right], u[i_result] );
|
||||
|
||||
// back propagate S(x) * f''(x) * R
|
||||
// (here is where we must check for cross terms)
|
||||
if( s[i_result] & vx[i_left] & vx[i_right] )
|
||||
{ v[i_left] = set_union(v[i_left], r[i_right] );
|
||||
v[i_right] = set_union(v[i_right], r[i_left] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
// bool sparsity
|
||||
virtual bool rev_sparse_hes(
|
||||
const vector<bool>& vx,
|
||||
const vector<bool>& s ,
|
||||
vector<bool>& t ,
|
||||
size_t q ,
|
||||
const vector<bool>& r ,
|
||||
const vector<bool>& u ,
|
||||
vector<bool>& v ,
|
||||
const vector<double>& x )
|
||||
{
|
||||
size_t nr_left = size_t( CppAD::Integer( x[0] ) );
|
||||
size_t n_middle = size_t( CppAD::Integer( x[1] ) );
|
||||
size_t nc_right = size_t( CppAD::Integer( x[2] ) );
|
||||
size_t nx = 3 + (nr_left + nc_right) * n_middle;
|
||||
# ifndef NDEBUG
|
||||
size_t ny = nr_left * nc_right;
|
||||
# endif
|
||||
assert( x.size() == nx );
|
||||
assert( vx.size() == nx );
|
||||
assert( t.size() == nx );
|
||||
assert( r.size() == nx * q );
|
||||
assert( v.size() == nx * q );
|
||||
assert( s.size() == ny );
|
||||
assert( u.size() == ny * q );
|
||||
size_t i, j, p;
|
||||
//
|
||||
// initilaize sparsity patterns as false
|
||||
for(j = 0; j < nx; j++)
|
||||
{ t[j] = false;
|
||||
for(p = 0; p < q; p++)
|
||||
v[j * q + p] = false;
|
||||
}
|
||||
size_t nk = 1;
|
||||
size_t k = 0;
|
||||
for(i = 0; i < nr_left; i++)
|
||||
{ for(j = 0; j < nc_right; j++)
|
||||
{ size_t i_result = result(
|
||||
i, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
for(size_t ell = 0; ell < n_middle; ell++)
|
||||
{ size_t i_left = left(
|
||||
i, ell, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
size_t i_right = right(
|
||||
ell, j, k, nk, nr_left, n_middle, nc_right
|
||||
);
|
||||
//
|
||||
// Compute sparsity for T(x) = S(x) * f'(x).
|
||||
// We so not need to use vx with f'(x) propagation.
|
||||
t[i_left] |= bool( s[i_result] );
|
||||
t[i_right] |= bool( s[i_result] );
|
||||
|
||||
// V(x) = f'(x)^T * U(x) + S(x) * f''(x) * R
|
||||
// U(x) = g''(y) * f'(x) * R
|
||||
// S(x) = g'(y)
|
||||
|
||||
// back propagate f'(x)^T * U(x)
|
||||
// (no need to use vx with f'(x) propogation)
|
||||
for(p = 0; p < q; p++)
|
||||
{ v[ i_left * q + p] |= bool( u[ i_result * q + p] );
|
||||
v[ i_right * q + p] |= bool( u[ i_result * q + p] );
|
||||
}
|
||||
|
||||
// back propagate S(x) * f''(x) * R
|
||||
// (here is where we must check for cross terms)
|
||||
if( s[i_result] & vx[i_left] & vx[i_right] )
|
||||
{ for(p = 0; p < q; p++)
|
||||
{ v[i_left * q + p] |= bool( r[i_right * q + p] );
|
||||
v[i_right * q + p] |= bool( r[i_left * q + p] );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return true;
|
||||
}
|
||||
/* %$$
|
||||
$head End Class Definition$$
|
||||
$srccode%cpp% */
|
||||
}; // End of mat_mul class
|
||||
} // End empty namespace
|
||||
/* %$$
|
||||
$comment end nospell$$
|
||||
$end
|
||||
*/
|
||||
|
||||
|
||||
# endif
|
||||
Reference in New Issue
Block a user