21 #ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
22 #define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
24 #include <opm/common/utility/platform_dependent/disable_warnings.h>
26 #include <Eigen/Eigen>
27 #include <Eigen/Sparse>
28 #include <opm/autodiff/fastSparseOperations.hpp>
30 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
32 #include <opm/autodiff/AutoDiffMatrix.hpp>
94 template <
typename Scalar>
99 typedef Eigen::Array<Scalar, Eigen::Dynamic, 1>
V;
132 const int num_elem = val.
size();
133 const int num_blocks = blocksizes.size();
135 jac.resize(num_blocks);
136 for (
int i = 0; i < num_blocks; ++i) {
137 jac[i] =
M(num_elem, blocksizes[i]);
153 const int num_elem = val.
size();
154 const int num_blocks = blocksizes.size();
156 jac.resize(num_blocks);
157 for (
int i = 0; i < num_blocks; ++i) {
158 jac[i] =
M(num_elem, blocksizes[i]);
161 assert(blocksizes[index] == num_elem);
176 return variable(index, std::move(value), blocksizes);
197 std::vector<M> jac_copy(jac);
198 return AutoDiffBlock(std::move(val_copy), std::move(jac_copy));
203 static std::vector<AutoDiffBlock>
variables(
const std::vector<V>& initial_values)
205 const int num_vars = initial_values.size();
206 std::vector<int> bpat;
207 for (
int v = 0; v < num_vars; ++v) {
208 bpat.push_back(initial_values[v].
size());
210 std::vector<AutoDiffBlock> vars;
211 for (
int v = 0; v < num_vars; ++v) {
212 vars.emplace_back(
variable(v, initial_values[v], bpat));
222 }
else if (!rhs.jac_.empty()) {
227 #pragma omp parallel for schedule(static)
228 for (
int block = 0; block < num_blocks; ++block) {
229 assert(jac_[block].rows() == rhs.jac_[block].rows());
230 assert(jac_[block].cols() == rhs.jac_[block].cols());
231 jac_[block] += rhs.jac_[block];
245 jac_.resize(num_blocks);
246 #pragma omp parallel for schedule(static)
247 for (
int block = 0; block < num_blocks; ++block) {
248 jac_[block] = rhs.jac_[block] * (-1.0);
250 }
else if (!rhs.jac_.empty()) {
255 #pragma omp parallel for schedule(static)
256 for (
int block = 0; block < num_blocks; ++block) {
257 assert(jac_[block].rows() == rhs.jac_[block].rows());
258 assert(jac_[block].cols() == rhs.jac_[block].cols());
259 jac_[block] -= rhs.jac_[block];
271 if (jac_.empty() && rhs.jac_.empty()) {
277 if (rhs.jac_.empty()) {
278 return *
this + rhs.val_;
280 std::vector<M> jac = jac_;
283 #pragma omp parallel for schedule(static)
284 for (
int block = 0; block < num_blocks; ++block) {
285 assert(jac[block].rows() == rhs.jac_[block].rows());
286 assert(jac[block].cols() == rhs.jac_[block].cols());
287 jac[block] += rhs.jac_[block];
289 return function(val_ + rhs.val_, std::move(jac));
295 if (jac_.empty() && rhs.jac_.empty()) {
301 if (rhs.jac_.empty()) {
302 return *
this - rhs.val_;
304 std::vector<M> jac = jac_;
307 #pragma omp parallel for schedule(static)
308 for (
int block = 0; block < num_blocks; ++block) {
309 assert(jac[block].rows() == rhs.jac_[block].rows());
310 assert(jac[block].cols() == rhs.jac_[block].cols());
311 jac[block] -= rhs.jac_[block];
313 return function(val_ - rhs.val_, std::move(jac));
319 if (jac_.empty() && rhs.jac_.empty()) {
325 if (rhs.jac_.empty()) {
326 return *
this * rhs.val_;
329 std::vector<M> jac(num_blocks);
331 M D1(val_.matrix().asDiagonal());
332 M D2(rhs.val_.matrix().asDiagonal());
333 #pragma omp parallel for schedule(dynamic)
334 for (
int block = 0; block < num_blocks; ++block) {
335 assert(jac_[block].rows() == rhs.jac_[block].rows());
336 assert(jac_[block].cols() == rhs.jac_[block].cols());
337 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
338 jac[block] =
M( D2.rows(), jac_[block].cols() );
340 else if( jac_[block].nonZeros() == 0 )
341 jac[block] = D1*rhs.jac_[block];
342 else if ( rhs.jac_[block].nonZeros() == 0 ) {
343 jac[block] = D2*jac_[block];
346 jac[block] = D2*jac_[block];
347 jac[block] += D1*rhs.jac_[block];
350 return function(val_ * rhs.val_, std::move(jac));
356 if (jac_.empty() && rhs.jac_.empty()) {
362 if (rhs.jac_.empty()) {
363 return *
this / rhs.val_;
366 std::vector<M> jac(num_blocks);
368 M D1(val_.matrix().asDiagonal());
369 M D2(rhs.val_.matrix().asDiagonal());
370 M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
371 #pragma omp parallel for schedule(dynamic)
372 for (
int block = 0; block < num_blocks; ++block) {
373 assert(jac_[block].rows() == rhs.jac_[block].rows());
374 assert(jac_[block].cols() == rhs.jac_[block].cols());
375 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
376 jac[block] =
M( D3.rows(), jac_[block].cols() );
378 else if( jac_[block].nonZeros() == 0 ) {
379 jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0);
381 else if ( rhs.jac_[block].nonZeros() == 0 )
383 jac[block] = D3 * (D2*jac_[block]);
386 jac[block] = D3 * (D2*jac_[block] + (D1*rhs.jac_[block]*(-1.0)));
389 return function(val_ / rhs.val_, std::move(jac));
393 template <
class Ostream>
397 int num_blocks = jac_.size();
398 os <<
"Value =\n" << val_ <<
"\n\nJacobian =\n";
399 for (
int i = 0; i < num_blocks; ++i) {
400 Eigen::SparseMatrix<double> m;
402 os <<
"Sub Jacobian #" << i <<
'\n' << m <<
"\n";
410 val_.swap(other.val_);
411 jac_.swap(other.jac_);
430 std::vector<int> bp(nb);
431 for (
int block = 0; block < nb; ++block) {
432 bp[block] = jac_[block].cols();
455 AutoDiffBlock(
V&& val)
456 : val_(std::move(val))
460 AutoDiffBlock(
V&& val, std::vector<M>&& jac)
461 : val_(std::move(val)), jac_(std::move(jac))
464 const int num_elem = val_.size();
465 const int num_blocks = jac_.size();
466 for (
int block = 0; block < num_blocks; ++block) {
467 assert(num_elem == jac_[block].rows());
480 template <
class Ostream,
typename Scalar>
482 operator<<(Ostream& os, const AutoDiffBlock<Scalar>& fw)
489 template <
typename Scalar>
494 std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
495 assert(lhs.
cols() == rhs.
value().rows());
496 #pragma omp parallel for schedule(dynamic)
497 for (
int block = 0; block < num_blocks; ++block) {
506 template <
typename Scalar>
511 std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
512 assert(lhs.cols() == rhs.
value().rows());
513 for (
int block = 0; block < num_blocks; ++block) {
522 template <
typename Scalar>
531 template <
typename Scalar>
540 template <
typename Scalar>
549 template <
typename Scalar>
558 template <
typename Scalar>
567 template <
typename Scalar>
576 template <
typename Scalar>
585 template <
typename Scalar>
600 template <
typename Scalar>
604 std::vector< typename AutoDiffBlock<Scalar>::M > jac;
606 for (
int block=0; block<lhs.
numBlocks(); block++) {
607 jac.emplace_back( lhs.
derivative()[block] * rhs );
620 template <
typename Scalar>
635 template <
typename Scalar>
637 const double exponent)
643 std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.
numBlocks());
644 for (
int block = 0; block < base.
numBlocks(); block++) {
658 template <
typename Scalar>
673 template <
typename Scalar>
688 template <
typename Scalar>
692 const int num_elem = base.
value().size();
693 assert(exponent.
size() == num_elem);
695 for (
int i = 0; i < num_elem; ++i) {
696 val[i] = std::pow(base.
value()[i], exponent.
value()[i]);
707 std::vector< typename AutoDiffBlock<Scalar>::M > jac (num_blocks);
711 for (
int i = 0; i < num_elem; ++i) {
712 der1[i] *= std::log(base.
value()[i]);
714 std::vector< typename AutoDiffBlock<Scalar>::M > jac1 (exponent.
numBlocks());
716 for (
int block = 0; block < exponent.
numBlocks(); block++) {
718 jac[block] = jac1[block];
724 for (
int i = 0; i < num_elem; ++i) {
725 der2[i] *= std::pow(base.
value()[i], exponent.
value()[i] - 1.0);
727 std::vector< typename AutoDiffBlock<Scalar>::M > jac2 (base.
numBlocks());
729 for (
int block = 0; block < base.
numBlocks(); block++) {
732 jac[block] += jac2[block];
734 jac[block] = jac2[block];
748 #endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED
void swap(AutoDiffBlock &other)
Efficient swap function.
Definition: AutoDiffBlock.hpp:408
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
const std::vector< M > & derivative() const
Function derivatives.
Definition: AutoDiffBlock.hpp:444
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
AutoDiffBlock operator/(const AutoDiffBlock &rhs) const
Elementwise operator /.
Definition: AutoDiffBlock.hpp:354
static AutoDiffMatrix createIdentity(const int num_rows_cols)
Creates an identity matrix with num_rows_cols x num_rows_cols entries.
Definition: AutoDiffMatrix.hpp:81
AutoDiffBlock & operator+=(const AutoDiffBlock &rhs)
Elementwise operator +=.
Definition: AutoDiffBlock.hpp:218
static AutoDiffBlock variable(const int index, const V &val, const std::vector< int > &blocksizes)
Create an AutoDiffBlock representing a single variable block.
Definition: AutoDiffBlock.hpp:173
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:101
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
A class for forward-mode automatic differentiation with vector values and sparse jacobian matrices...
Definition: AutoDiffBlock.hpp:95
int numBlocks() const
Number of Jacobian blocks.
Definition: AutoDiffBlock.hpp:421
AutoDiffBlock operator*(const AutoDiffBlock &rhs) const
Elementwise operator *.
Definition: AutoDiffBlock.hpp:317
AutoDiffBlock< Scalar > pow(const AutoDiffBlock< Scalar > &base, const double exponent)
Computes the value of base raised to the power of exponent.
Definition: AutoDiffBlock.hpp:636
int cols() const
Returns number of columns in the matrix.
Definition: AutoDiffMatrix.hpp:575
Ostream & print(Ostream &os) const
I/O.
Definition: AutoDiffBlock.hpp:395
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
AutoDiffBlock operator-(const AutoDiffBlock &rhs) const
Elementwise operator -.
Definition: AutoDiffBlock.hpp:293
std::vector< int > blockPattern() const
Sizes (number of columns) of Jacobian blocks.
Definition: AutoDiffBlock.hpp:427
static AutoDiffBlock constant(const V &val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:118
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
static AutoDiffBlock variable(const int index, V &&val, const std::vector< int > &blocksizes)
Create an AutoDiffBlock representing a single variable block.
Definition: AutoDiffBlock.hpp:150
AutoDiffBlock operator+(const AutoDiffBlock &rhs) const
Elementwise operator +.
Definition: AutoDiffBlock.hpp:269
static AutoDiffBlock constant(const V &val, const std::vector< int > &blocksizes)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:129
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Construct a set of primary variables, each initialized to a given vector.
Definition: AutoDiffBlock.hpp:203
AutoDiffBlock & operator-=(const AutoDiffBlock &rhs)
Elementwise operator -=.
Definition: AutoDiffBlock.hpp:241