20 #ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
21 #define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
23 #include <opm/common/utility/platform_dependent/disable_warnings.h>
25 #include <Eigen/Eigen>
26 #include <Eigen/Sparse>
29 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
93 template <
typename Scalar>
98 typedef Eigen::Array<Scalar, Eigen::Dynamic, 1>
V;
131 const int num_elem = val.
size();
132 const int num_blocks = blocksizes.size();
134 jac.resize(num_blocks);
135 for (
int i = 0; i < num_blocks; ++i) {
136 jac[i] =
M(num_elem, blocksizes[i]);
152 const int num_elem = val.
size();
153 const int num_blocks = blocksizes.size();
155 jac.resize(num_blocks);
156 for (
int i = 0; i < num_blocks; ++i) {
157 jac[i] =
M(num_elem, blocksizes[i]);
160 assert(blocksizes[index] == num_elem);
175 return variable(index, std::move(value), blocksizes);
196 std::vector<M> jac_copy(jac);
197 return AutoDiffBlock(std::move(val_copy), std::move(jac_copy));
202 static std::vector<AutoDiffBlock>
variables(
const std::vector<V>& initial_values)
204 const int num_vars = initial_values.size();
205 std::vector<int> bpat;
206 for (
int v = 0; v < num_vars; ++v) {
207 bpat.push_back(initial_values[v].
size());
209 std::vector<AutoDiffBlock> vars;
210 for (
int v = 0; v < num_vars; ++v) {
211 vars.emplace_back(
variable(v, initial_values[v], bpat));
221 }
else if (!rhs.jac_.empty()) {
226 #pragma omp parallel for schedule(static)
227 for (
int block = 0; block < num_blocks; ++block) {
228 assert(jac_[block].rows() == rhs.jac_[block].rows());
229 assert(jac_[block].cols() == rhs.jac_[block].cols());
230 jac_[block] += rhs.jac_[block];
244 jac_.resize(num_blocks);
245 #pragma omp parallel for schedule(static)
246 for (
int block = 0; block < num_blocks; ++block) {
247 jac_[block] = rhs.jac_[block] * (-1.0);
249 }
else if (!rhs.jac_.empty()) {
254 #pragma omp parallel for schedule(static)
255 for (
int block = 0; block < num_blocks; ++block) {
256 assert(jac_[block].rows() == rhs.jac_[block].rows());
257 assert(jac_[block].cols() == rhs.jac_[block].cols());
258 jac_[block] -= rhs.jac_[block];
270 if (jac_.empty() && rhs.jac_.empty()) {
276 if (rhs.jac_.empty()) {
277 return *
this + rhs.val_;
279 std::vector<M> jac = jac_;
282 #pragma omp parallel for schedule(static)
283 for (
int block = 0; block < num_blocks; ++block) {
284 assert(jac[block].rows() == rhs.jac_[block].rows());
285 assert(jac[block].cols() == rhs.jac_[block].cols());
286 jac[block] += rhs.jac_[block];
288 return function(val_ + rhs.val_, std::move(jac));
294 if (jac_.empty() && rhs.jac_.empty()) {
300 if (rhs.jac_.empty()) {
301 return *
this - rhs.val_;
303 std::vector<M> jac = jac_;
306 #pragma omp parallel for schedule(static)
307 for (
int block = 0; block < num_blocks; ++block) {
308 assert(jac[block].rows() == rhs.jac_[block].rows());
309 assert(jac[block].cols() == rhs.jac_[block].cols());
310 jac[block] -= rhs.jac_[block];
312 return function(val_ - rhs.val_, std::move(jac));
318 if (jac_.empty() && rhs.jac_.empty()) {
324 if (rhs.jac_.empty()) {
325 return *
this * rhs.val_;
328 std::vector<M> jac(num_blocks);
330 M D1(val_.matrix().asDiagonal());
331 M D2(rhs.val_.matrix().asDiagonal());
332 #pragma omp parallel for schedule(dynamic)
333 for (
int block = 0; block < num_blocks; ++block) {
334 assert(jac_[block].rows() == rhs.jac_[block].rows());
335 assert(jac_[block].cols() == rhs.jac_[block].cols());
336 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
337 jac[block] =
M( D2.rows(), jac_[block].cols() );
339 else if( jac_[block].nonZeros() == 0 )
340 jac[block] = D1*rhs.jac_[block];
341 else if ( rhs.jac_[block].nonZeros() == 0 ) {
342 jac[block] = D2*jac_[block];
345 jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
348 return function(val_ * rhs.val_, std::move(jac));
354 if (jac_.empty() && rhs.jac_.empty()) {
360 if (rhs.jac_.empty()) {
361 return *
this / rhs.val_;
364 std::vector<M> jac(num_blocks);
366 M D1(val_.matrix().asDiagonal());
367 M D2(rhs.val_.matrix().asDiagonal());
368 M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
369 #pragma omp parallel for schedule(dynamic)
370 for (
int block = 0; block < num_blocks; ++block) {
371 assert(jac_[block].rows() == rhs.jac_[block].rows());
372 assert(jac_[block].cols() == rhs.jac_[block].cols());
373 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
374 jac[block] =
M( D3.rows(), jac_[block].cols() );
376 else if( jac_[block].nonZeros() == 0 ) {
377 jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0);
379 else if ( rhs.jac_[block].nonZeros() == 0 )
381 jac[block] = D3 * (D2*jac_[block]);
384 jac[block] = D3 * (D2*jac_[block] + (D1*rhs.jac_[block]*(-1.0)));
387 return function(val_ / rhs.val_, std::move(jac));
391 template <
class Ostream>
395 int num_blocks = jac_.size();
396 os <<
"Value =\n" << val_ <<
"\n\nJacobian =\n";
397 for (
int i = 0; i < num_blocks; ++i) {
398 Eigen::SparseMatrix<double> m;
400 os <<
"Sub Jacobian #" << i <<
'\n' << m <<
"\n";
408 val_.swap(other.val_);
409 jac_.swap(other.jac_);
428 std::vector<int> bp(nb);
429 for (
int block = 0; block < nb; ++block) {
430 bp[block] = jac_[block].cols();
453 AutoDiffBlock(V&& val)
454 : val_(
std::move(val))
458 AutoDiffBlock(V&& val, std::vector<M>&& jac)
459 : val_(
std::move(val)), jac_(
std::move(jac))
462 const int num_elem = val_.size();
463 const int num_blocks = jac_.size();
464 for (
int block = 0; block < num_blocks; ++block) {
465 assert(num_elem == jac_[block].rows());
478 template <
class Ostream,
typename Scalar>
480 operator<<(Ostream& os, const AutoDiffBlock<Scalar>& fw)
487 template <
typename Scalar>
493 assert(lhs.
cols() == rhs.
value().rows());
494 #pragma omp parallel for schedule(dynamic)
495 for (
int block = 0; block < num_blocks; ++block) {
504 template <
typename Scalar>
510 assert(lhs.cols() == rhs.
value().rows());
511 for (
int block = 0; block < num_blocks; ++block) {
520 template <
typename Scalar>
529 template <
typename Scalar>
538 template <
typename Scalar>
547 template <
typename Scalar>
556 template <
typename Scalar>
565 template <
typename Scalar>
574 template <
typename Scalar>
583 template <
typename Scalar>
598 template <
typename Scalar>
604 for (
int block=0; block<lhs.
numBlocks(); block++) {
605 jac.emplace_back( lhs.
derivative()[block] * rhs );
618 template <
typename Scalar>
630 #endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED
AutoDiff< Scalar > operator+(const AutoDiff< Scalar > &lhs, const AutoDiff< Scalar > &rhs)
Definition: AutoDiff.hpp:148
AutoDiff< Scalar > operator*(const AutoDiff< Scalar > &lhs, const AutoDiff< Scalar > &rhs)
Definition: AutoDiff.hpp:211
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:98
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Definition: AutoDiffBlock.hpp:202
static AutoDiffMatrix createIdentity(const int num_rows_cols)
Definition: AutoDiffMatrix.hpp:81
Ostream & print(Ostream &os) const
I/O.
Definition: AutoDiffBlock.hpp:393
Definition: AdditionalObjectDeleter.hpp:22
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:100
std::vector< int > blockPattern() const
Sizes (number of columns) of Jacobian blocks.
Definition: AutoDiffBlock.hpp:425
int cols() const
Definition: AutoDiffMatrix.hpp:563
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:103
static AutoDiffBlock constant(const V &val)
Definition: AutoDiffBlock.hpp:117
int numBlocks() const
Number of Jacobian blocks.
Definition: AutoDiffBlock.hpp:419
AutoDiffBlock operator/(const AutoDiffBlock &rhs) const
Elementwise operator /.
Definition: AutoDiffBlock.hpp:352
AutoDiffBlock & operator-=(const AutoDiffBlock &rhs)
Elementwise operator -=.
Definition: AutoDiffBlock.hpp:240
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Definition: AutoDiffMatrix.hpp:696
AutoDiffBlock operator*(const AutoDiffBlock &rhs) const
Elementwise operator *.
Definition: AutoDiffBlock.hpp:316
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:413
AutoDiffBlock & operator+=(const AutoDiffBlock &rhs)
Elementwise operator +=.
Definition: AutoDiffBlock.hpp:217
AutoDiffBlock operator+(const AutoDiffBlock &rhs) const
Elementwise operator +.
Definition: AutoDiffBlock.hpp:268
static AutoDiffBlock constant(V &&val)
Definition: AutoDiffBlock.hpp:110
AutoDiff< Scalar > operator/(const AutoDiff< Scalar > &lhs, const AutoDiff< Scalar > &rhs)
Definition: AutoDiff.hpp:244
Definition: AutoDiffMatrix.hpp:43
static AutoDiffBlock constant(const V &val, const std::vector< int > &blocksizes)
Definition: AutoDiffBlock.hpp:128
AutoDiff< Scalar > operator-(const AutoDiff< Scalar > &lhs, const AutoDiff< Scalar > &rhs)
Definition: AutoDiff.hpp:181
void swap(AutoDiffBlock &other)
Efficient swap function.
Definition: AutoDiffBlock.hpp:406
static AutoDiffBlock variable(const int index, const V &val, const std::vector< int > &blocksizes)
Definition: AutoDiffBlock.hpp:172
ADB::M M
Definition: BlackoilModelBase_impl.hpp:85
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Definition: AutoDiffBlock.hpp:183
AutoDiffBlock operator-(const AutoDiffBlock &rhs) const
Elementwise operator -.
Definition: AutoDiffBlock.hpp:292
Definition: AutoDiffBlock.hpp:94
static AutoDiffBlock variable(const int index, V &&val, const std::vector< int > &blocksizes)
Definition: AutoDiffBlock.hpp:149
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:436
const std::vector< M > & derivative() const
Function derivatives.
Definition: AutoDiffBlock.hpp:442