Added ForwardBlock operator* for scalars.

This commit is contained in:
Jens Olav Nygaard 2013-08-02 12:24:12 +02:00
parent 23e2034118
commit d3a02e4891
2 changed files with 90 additions and 33 deletions

View File

@ -242,7 +242,17 @@ namespace AutoDiff
V val_;
std::vector<M> jac_;
};
template <typename Sclr>
friend
ForwardBlock<Sclr> operator*(const ForwardBlock<Sclr> &lhs,
const Sclr &rhs);
template <typename Sclr>
friend
ForwardBlock<Sclr> operator*(const Sclr &lhs,
const ForwardBlock<Sclr> &rhs);
};
template <class Ostream, typename Scalar>
@ -333,6 +343,41 @@ namespace AutoDiff
}
/**
* @brief Operator for multiplication with a scalar on the right-hand side
*
* @param lhs The left-hand side AD forward block
* @param rhs The scalar to multiply with
* @return The product
*/
template <typename Scalar>
ForwardBlock<Scalar> operator*(const ForwardBlock<Scalar> &lhs,
const Scalar &rhs)
{
std::vector< Eigen::SparseMatrix<Scalar> > jac = lhs.jac_;
for (int block=0; block<lhs.numBlocks(); block++) {
jac[block] = lhs.jac_[block] * rhs;
}
auto val = lhs.val_ * rhs;
return ForwardBlock<Scalar>::function(val, jac);
}
/**
* @brief Operator for multiplication with a scalar on the left-hand side
*
* @param lhs The scalar to multiply with
* @param rhs The right-hand side AD forward block
* @return The product
*/
template <typename Scalar>
ForwardBlock<Scalar> operator*(const Scalar &lhs,
const ForwardBlock<Scalar> &rhs)
{
return rhs * lhs; // Commutative operation.
}
} // namespace Autodiff

View File

@ -1,53 +1,65 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE ScalarMultTest
#define BOOST_TEST_MODULE AutoDiffBlockTest
#include <opm/autodiff/AutoDiff.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <cmath>
#include <boost/test/unit_test.hpp>
#include <iostream>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
BOOST_AUTO_TEST_CASE(ScalarMultiplication)
{
typedef AutoDiff::ForwardBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
std::cout << "er her" << std::endl;
ADB::V vx(3);
vx << 1.0, 2.0, 10.0;
enum { FirstVar = 0, SecondVar = 1, ThirdVar = 2 };
typedef AutoDiff::Forward<double> AdFW;
ADB x = ADB::variable(FirstVar, vx, blocksizes);
const double atol = 1.0e-14;
ADB::V const_vector(3);
const_vector << 3.14, 3.14, 3.14;
// std::cout << "const_vector:\n" << const_vector << std::endl;
const ADB x2 = x * const_vector;
const ADB x3 = const_vector * x;
BOOST_CHECK_EQUAL( x2.value().cwiseNotEqual( x3.value() ).count(), 0 );
AdFW a = AdFW::variable(0.0);
AdFW b = AdFW::variable(1.0);
// The new operator:
AdFW two_a = a + a;
BOOST_CHECK_CLOSE(two_a.val(), 2*a.val(), atol);
BOOST_CHECK_CLOSE(two_a.der(), 2*a.der(), atol);
const ADB y2 = x * 3.14;
BOOST_CHECK_EQUAL( x2.value().cwiseNotEqual( y2.value() ).count(), 0 );
double av = a.val();
double ad = a.der();
a += b;
BOOST_CHECK_CLOSE(a.val(), av + b.val(), atol);
BOOST_CHECK_CLOSE(a.der(), ad + b.der(), atol);
av = a.val();
ad = a.der();
a += 1;
BOOST_CHECK_CLOSE(a.val(), av + 1, atol);
BOOST_CHECK_CLOSE(a.der(), ad , atol);
AdFW bpo = b + 1; // b plus one
BOOST_CHECK_CLOSE(bpo.val(), b.val() + 1, atol);
BOOST_CHECK_CLOSE(bpo.der(), b.der() , atol);
AdFW opb = 1 + b; // one plus b
BOOST_CHECK_CLOSE(opb.val(), b.val() + 1, atol);
BOOST_CHECK_CLOSE(opb.der(), b.der() , atol);
const ADB y3 = 3.14 * x;
BOOST_CHECK_EQUAL( x3.value().cwiseNotEqual( y3.value() ).count(), 0 );
}