00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
00022 #define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
00023
00024 #include <opm/common/utility/platform_dependent/disable_warnings.h>
00025
00026 #include <Eigen/Eigen>
00027 #include <Eigen/Sparse>
00028 #include <opm/autodiff/fastSparseOperations.hpp>
00029
00030 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
00031
00032 #include <opm/autodiff/AutoDiffMatrix.hpp>
00033
00034
00035 #include <utility>
00036 #include <vector>
00037 #include <cassert>
00038 #include <iostream>
00039
00040 namespace Opm
00041 {
00042
00094 template <typename Scalar>
00095 class AutoDiffBlock
00096 {
00097 public:
00099 typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
00101 typedef AutoDiffMatrix M;
00102
00104 static AutoDiffBlock null()
00105 {
00106 return AutoDiffBlock(V(), {});
00107 }
00108
00111 static AutoDiffBlock constant(V&& val)
00112 {
00113 return AutoDiffBlock(std::move(val));
00114 }
00115
00118 static AutoDiffBlock constant(const V& val)
00119 {
00120 return AutoDiffBlock(val);
00121 }
00122
00129 static AutoDiffBlock constant(const V& val, const std::vector<int>& blocksizes)
00130 {
00131 std::vector<M> jac;
00132 const int num_elem = val.size();
00133 const int num_blocks = blocksizes.size();
00134
00135 jac.resize(num_blocks);
00136 for (int i = 0; i < num_blocks; ++i) {
00137 jac[i] = M(num_elem, blocksizes[i]);
00138 }
00139 V val_copy(val);
00140 return AutoDiffBlock(std::move(val_copy), std::move(jac));
00141 }
00142
00150 static AutoDiffBlock variable(const int index, V&& val, const std::vector<int>& blocksizes)
00151 {
00152 std::vector<M> jac;
00153 const int num_elem = val.size();
00154 const int num_blocks = blocksizes.size();
00155
00156 jac.resize(num_blocks);
00157 for (int i = 0; i < num_blocks; ++i) {
00158 jac[i] = M(num_elem, blocksizes[i]);
00159 }
00160
00161 assert(blocksizes[index] == num_elem);
00162 jac[index] = M::createIdentity(val.size());
00163 return AutoDiffBlock(std::move(val), std::move(jac));
00164 }
00165
00173 static AutoDiffBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
00174 {
00175 V value = val;
00176 return variable(index, std::move(value), blocksizes);
00177 }
00178
00184 static AutoDiffBlock function(V&& val, std::vector<M>&& jac)
00185 {
00186 return AutoDiffBlock(std::move(val), std::move(jac));
00187 }
00188
00194 static AutoDiffBlock function(const V& val, const std::vector<M>& jac)
00195 {
00196 V val_copy(val);
00197 std::vector<M> jac_copy(jac);
00198 return AutoDiffBlock(std::move(val_copy), std::move(jac_copy));
00199 }
00200
00203 static std::vector<AutoDiffBlock> variables(const std::vector<V>& initial_values)
00204 {
00205 const int num_vars = initial_values.size();
00206 std::vector<int> bpat;
00207 for (int v = 0; v < num_vars; ++v) {
00208 bpat.push_back(initial_values[v].size());
00209 }
00210 std::vector<AutoDiffBlock> vars;
00211 for (int v = 0; v < num_vars; ++v) {
00212 vars.emplace_back(variable(v, initial_values[v], bpat));
00213 }
00214 return vars;
00215 }
00216
00218 AutoDiffBlock& operator+=(const AutoDiffBlock& rhs)
00219 {
00220 if (jac_.empty()) {
00221 jac_ = rhs.jac_;
00222 } else if (!rhs.jac_.empty()) {
00223 assert (numBlocks() == rhs.numBlocks());
00224 assert (value().size() == rhs.value().size());
00225
00226 const int num_blocks = numBlocks();
00227 #pragma omp parallel for schedule(static)
00228 for (int block = 0; block < num_blocks; ++block) {
00229 assert(jac_[block].rows() == rhs.jac_[block].rows());
00230 assert(jac_[block].cols() == rhs.jac_[block].cols());
00231 jac_[block] += rhs.jac_[block];
00232 }
00233 }
00234
00235 val_ += rhs.val_;
00236
00237 return *this;
00238 }
00239
00241 AutoDiffBlock& operator-=(const AutoDiffBlock& rhs)
00242 {
00243 if (jac_.empty()) {
00244 const int num_blocks = rhs.numBlocks();
00245 jac_.resize(num_blocks);
00246 #pragma omp parallel for schedule(static)
00247 for (int block = 0; block < num_blocks; ++block) {
00248 jac_[block] = rhs.jac_[block] * (-1.0);
00249 }
00250 } else if (!rhs.jac_.empty()) {
00251 assert (numBlocks() == rhs.numBlocks());
00252 assert (value().size() == rhs.value().size());
00253
00254 const int num_blocks = numBlocks();
00255 #pragma omp parallel for schedule(static)
00256 for (int block = 0; block < num_blocks; ++block) {
00257 assert(jac_[block].rows() == rhs.jac_[block].rows());
00258 assert(jac_[block].cols() == rhs.jac_[block].cols());
00259 jac_[block] -= rhs.jac_[block];
00260 }
00261 }
00262
00263 val_ -= rhs.val_;
00264
00265 return *this;
00266 }
00267
00269 AutoDiffBlock operator+(const AutoDiffBlock& rhs) const
00270 {
00271 if (jac_.empty() && rhs.jac_.empty()) {
00272 return constant(val_ + rhs.val_);
00273 }
00274 if (jac_.empty()) {
00275 return val_ + rhs;
00276 }
00277 if (rhs.jac_.empty()) {
00278 return *this + rhs.val_;
00279 }
00280 std::vector<M> jac = jac_;
00281 assert(numBlocks() == rhs.numBlocks());
00282 int num_blocks = numBlocks();
00283 #pragma omp parallel for schedule(static)
00284 for (int block = 0; block < num_blocks; ++block) {
00285 assert(jac[block].rows() == rhs.jac_[block].rows());
00286 assert(jac[block].cols() == rhs.jac_[block].cols());
00287 jac[block] += rhs.jac_[block];
00288 }
00289 return function(val_ + rhs.val_, std::move(jac));
00290 }
00291
00293 AutoDiffBlock operator-(const AutoDiffBlock& rhs) const
00294 {
00295 if (jac_.empty() && rhs.jac_.empty()) {
00296 return constant(val_ - rhs.val_);
00297 }
00298 if (jac_.empty()) {
00299 return val_ - rhs;
00300 }
00301 if (rhs.jac_.empty()) {
00302 return *this - rhs.val_;
00303 }
00304 std::vector<M> jac = jac_;
00305 assert(numBlocks() == rhs.numBlocks());
00306 int num_blocks = numBlocks();
00307 #pragma omp parallel for schedule(static)
00308 for (int block = 0; block < num_blocks; ++block) {
00309 assert(jac[block].rows() == rhs.jac_[block].rows());
00310 assert(jac[block].cols() == rhs.jac_[block].cols());
00311 jac[block] -= rhs.jac_[block];
00312 }
00313 return function(val_ - rhs.val_, std::move(jac));
00314 }
00315
00317 AutoDiffBlock operator*(const AutoDiffBlock& rhs) const
00318 {
00319 if (jac_.empty() && rhs.jac_.empty()) {
00320 return constant(val_ * rhs.val_);
00321 }
00322 if (jac_.empty()) {
00323 return val_ * rhs;
00324 }
00325 if (rhs.jac_.empty()) {
00326 return *this * rhs.val_;
00327 }
00328 int num_blocks = numBlocks();
00329 std::vector<M> jac(num_blocks);
00330 assert(numBlocks() == rhs.numBlocks());
00331 M D1(val_.matrix().asDiagonal());
00332 M D2(rhs.val_.matrix().asDiagonal());
00333 #pragma omp parallel for schedule(dynamic)
00334 for (int block = 0; block < num_blocks; ++block) {
00335 assert(jac_[block].rows() == rhs.jac_[block].rows());
00336 assert(jac_[block].cols() == rhs.jac_[block].cols());
00337 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
00338 jac[block] = M( D2.rows(), jac_[block].cols() );
00339 }
00340 else if( jac_[block].nonZeros() == 0 )
00341 jac[block] = D1*rhs.jac_[block];
00342 else if ( rhs.jac_[block].nonZeros() == 0 ) {
00343 jac[block] = D2*jac_[block];
00344 }
00345 else {
00346 jac[block] = D2*jac_[block];
00347 jac[block] += D1*rhs.jac_[block];
00348 }
00349 }
00350 return function(val_ * rhs.val_, std::move(jac));
00351 }
00352
00354 AutoDiffBlock operator/(const AutoDiffBlock& rhs) const
00355 {
00356 if (jac_.empty() && rhs.jac_.empty()) {
00357 return constant(val_ / rhs.val_);
00358 }
00359 if (jac_.empty()) {
00360 return val_ / rhs;
00361 }
00362 if (rhs.jac_.empty()) {
00363 return *this / rhs.val_;
00364 }
00365 int num_blocks = numBlocks();
00366 std::vector<M> jac(num_blocks);
00367 assert(numBlocks() == rhs.numBlocks());
00368 M D1(val_.matrix().asDiagonal());
00369 M D2(rhs.val_.matrix().asDiagonal());
00370 M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
00371 #pragma omp parallel for schedule(dynamic)
00372 for (int block = 0; block < num_blocks; ++block) {
00373 assert(jac_[block].rows() == rhs.jac_[block].rows());
00374 assert(jac_[block].cols() == rhs.jac_[block].cols());
00375 if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
00376 jac[block] = M( D3.rows(), jac_[block].cols() );
00377 }
00378 else if( jac_[block].nonZeros() == 0 ) {
00379 jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0);
00380 }
00381 else if ( rhs.jac_[block].nonZeros() == 0 )
00382 {
00383 jac[block] = D3 * (D2*jac_[block]);
00384 }
00385 else {
00386 jac[block] = D3 * (D2*jac_[block] + (D1*rhs.jac_[block]*(-1.0)));
00387 }
00388 }
00389 return function(val_ / rhs.val_, std::move(jac));
00390 }
00391
00393 template <class Ostream>
00394 Ostream&
00395 print(Ostream& os) const
00396 {
00397 int num_blocks = jac_.size();
00398 os << "Value =\n" << val_ << "\n\nJacobian =\n";
00399 for (int i = 0; i < num_blocks; ++i) {
00400 Eigen::SparseMatrix<double> m;
00401 jac_[i].toSparse(m);
00402 os << "Sub Jacobian #" << i << '\n' << m << "\n";
00403 }
00404 return os;
00405 }
00406
00408 void swap(AutoDiffBlock& other)
00409 {
00410 val_.swap(other.val_);
00411 jac_.swap(other.jac_);
00412 }
00413
00415 int size() const
00416 {
00417 return val_.size();
00418 }
00419
00421 int numBlocks() const
00422 {
00423 return jac_.size();
00424 }
00425
00427 std::vector<int> blockPattern() const
00428 {
00429 const int nb = numBlocks();
00430 std::vector<int> bp(nb);
00431 for (int block = 0; block < nb; ++block) {
00432 bp[block] = jac_[block].cols();
00433 }
00434 return bp;
00435 }
00436
00438 const V& value() const
00439 {
00440 return val_;
00441 }
00442
00444 const std::vector<M>& derivative() const
00445 {
00446 return jac_;
00447 }
00448
00449 private:
00450 AutoDiffBlock(const V& val)
00451 : val_(val)
00452 {
00453 }
00454
00455 AutoDiffBlock(V&& val)
00456 : val_(std::move(val))
00457 {
00458 }
00459
00460 AutoDiffBlock(V&& val, std::vector<M>&& jac)
00461 : val_(std::move(val)), jac_(std::move(jac))
00462 {
00463 #ifndef NDEBUG
00464 const int num_elem = val_.size();
00465 const int num_blocks = jac_.size();
00466 for (int block = 0; block < num_blocks; ++block) {
00467 assert(num_elem == jac_[block].rows());
00468 }
00469 #endif
00470 }
00471
00472 V val_;
00473 std::vector<M> jac_;
00474 };
00475
00476
00477
00478
00480 template <class Ostream, typename Scalar>
00481 Ostream&
00482 operator<<(Ostream& os, const AutoDiffBlock<Scalar>& fw)
00483 {
00484 return fw.print(os);
00485 }
00486
00487
00489 template <typename Scalar>
00490 AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::M& lhs,
00491 const AutoDiffBlock<Scalar>& rhs)
00492 {
00493 int num_blocks = rhs.numBlocks();
00494 std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
00495 assert(lhs.cols() == rhs.value().rows());
00496 #pragma omp parallel for schedule(dynamic)
00497 for (int block = 0; block < num_blocks; ++block) {
00498 fastSparseProduct(lhs, rhs.derivative()[block], jac[block]);
00499 }
00500 typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix();
00501 return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
00502 }
00503
00504
00506 template <typename Scalar>
00507 AutoDiffBlock<Scalar> operator*(const Eigen::SparseMatrix<Scalar>& lhs,
00508 const AutoDiffBlock<Scalar>& rhs)
00509 {
00510 int num_blocks = rhs.numBlocks();
00511 std::vector<typename AutoDiffBlock<Scalar>::M> jac(num_blocks);
00512 assert(lhs.cols() == rhs.value().rows());
00513 for (int block = 0; block < num_blocks; ++block) {
00514 fastSparseProduct(lhs, rhs.derivative()[block], jac[block]);
00515 }
00516 typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix();
00517 return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
00518 }
00519
00520
00522 template <typename Scalar>
00523 AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::V& lhs,
00524 const AutoDiffBlock<Scalar>& rhs)
00525 {
00526 return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) * rhs;
00527 }
00528
00529
00531 template <typename Scalar>
00532 AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
00533 const typename AutoDiffBlock<Scalar>::V& rhs)
00534 {
00535 return rhs * lhs;
00536 }
00537
00538
00540 template <typename Scalar>
00541 AutoDiffBlock<Scalar> operator+(const typename AutoDiffBlock<Scalar>::V& lhs,
00542 const AutoDiffBlock<Scalar>& rhs)
00543 {
00544 return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) + rhs;
00545 }
00546
00547
00549 template <typename Scalar>
00550 AutoDiffBlock<Scalar> operator+(const AutoDiffBlock<Scalar>& lhs,
00551 const typename AutoDiffBlock<Scalar>::V& rhs)
00552 {
00553 return rhs + lhs;
00554 }
00555
00556
00558 template <typename Scalar>
00559 AutoDiffBlock<Scalar> operator-(const typename AutoDiffBlock<Scalar>::V& lhs,
00560 const AutoDiffBlock<Scalar>& rhs)
00561 {
00562 return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) - rhs;
00563 }
00564
00565
00567 template <typename Scalar>
00568 AutoDiffBlock<Scalar> operator-(const AutoDiffBlock<Scalar>& lhs,
00569 const typename AutoDiffBlock<Scalar>::V& rhs)
00570 {
00571 return lhs - AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
00572 }
00573
00574
00576 template <typename Scalar>
00577 AutoDiffBlock<Scalar> operator/(const typename AutoDiffBlock<Scalar>::V& lhs,
00578 const AutoDiffBlock<Scalar>& rhs)
00579 {
00580 return AutoDiffBlock<Scalar>::constant(lhs, rhs.blockPattern()) / rhs;
00581 }
00582
00583
00585 template <typename Scalar>
00586 AutoDiffBlock<Scalar> operator/(const AutoDiffBlock<Scalar>& lhs,
00587 const typename AutoDiffBlock<Scalar>::V& rhs)
00588 {
00589 return lhs / AutoDiffBlock<Scalar>::constant(rhs, lhs.blockPattern());
00590 }
00591
00592
00600 template <typename Scalar>
00601 AutoDiffBlock<Scalar> operator*(const AutoDiffBlock<Scalar>& lhs,
00602 const Scalar& rhs)
00603 {
00604 std::vector< typename AutoDiffBlock<Scalar>::M > jac;
00605 jac.reserve( lhs.numBlocks() );
00606 for (int block=0; block<lhs.numBlocks(); block++) {
00607 jac.emplace_back( lhs.derivative()[block] * rhs );
00608 }
00609 return AutoDiffBlock<Scalar>::function( lhs.value() * rhs, std::move(jac) );
00610 }
00611
00612
00620 template <typename Scalar>
00621 AutoDiffBlock<Scalar> operator*(const Scalar& lhs,
00622 const AutoDiffBlock<Scalar>& rhs)
00623 {
00624 return rhs * lhs;
00625 }
00626
00627
00635 template <typename Scalar>
00636 AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
00637 const double exponent)
00638 {
00639 const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exponent);
00640 const typename AutoDiffBlock<Scalar>::V derivative = exponent * base.value().pow(exponent - 1.0);
00641 const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
00642
00643 std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
00644 for (int block = 0; block < base.numBlocks(); block++) {
00645 fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]);
00646 }
00647
00648 return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
00649 }
00650
00658 template <typename Scalar>
00659 AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
00660 const typename AutoDiffBlock<Scalar>::V& exponent)
00661 {
00662
00663 return pow(base,AutoDiffBlock<Scalar>::constant(exponent));
00664 }
00665
00673 template <typename Scalar>
00674 AutoDiffBlock<Scalar> pow(const typename AutoDiffBlock<Scalar>::V& base,
00675 const AutoDiffBlock<Scalar>& exponent)
00676 {
00677
00678 return pow(AutoDiffBlock<Scalar>::constant(base),exponent);
00679 }
00680
00688 template <typename Scalar>
00689 AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
00690 const AutoDiffBlock<Scalar>& exponent)
00691 {
00692 const int num_elem = base.value().size();
00693 assert(exponent.size() == num_elem);
00694 typename AutoDiffBlock<Scalar>::V val (num_elem);
00695 for (int i = 0; i < num_elem; ++i) {
00696 val[i] = std::pow(base.value()[i], exponent.value()[i]);
00697 }
00698
00699
00700
00701
00702
00703 int num_blocks = std::max (base.numBlocks(), exponent.numBlocks());
00704 if (!base.derivative().empty() && !exponent.derivative().empty()) {
00705 assert(exponent.numBlocks() == base.numBlocks());
00706 }
00707 std::vector< typename AutoDiffBlock<Scalar>::M > jac (num_blocks);
00708
00709 if ( !exponent.derivative().empty() ) {
00710 typename AutoDiffBlock<Scalar>::V der1 = val;
00711 for (int i = 0; i < num_elem; ++i) {
00712 der1[i] *= std::log(base.value()[i]);
00713 }
00714 std::vector< typename AutoDiffBlock<Scalar>::M > jac1 (exponent.numBlocks());
00715 const typename AutoDiffBlock<Scalar>::M der1_diag(der1.matrix().asDiagonal());
00716 for (int block = 0; block < exponent.numBlocks(); block++) {
00717 fastSparseProduct(der1_diag, exponent.derivative()[block], jac1[block]);
00718 jac[block] = jac1[block];
00719 }
00720 }
00721
00722 if ( !base.derivative().empty() ) {
00723 typename AutoDiffBlock<Scalar>::V der2 = exponent.value();
00724 for (int i = 0; i < num_elem; ++i) {
00725 der2[i] *= std::pow(base.value()[i], exponent.value()[i] - 1.0);
00726 }
00727 std::vector< typename AutoDiffBlock<Scalar>::M > jac2 (base.numBlocks());
00728 const typename AutoDiffBlock<Scalar>::M der2_diag(der2.matrix().asDiagonal());
00729 for (int block = 0; block < base.numBlocks(); block++) {
00730 fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]);
00731 if (!exponent.derivative().empty()) {
00732 jac[block] += jac2[block];
00733 } else {
00734 jac[block] = jac2[block];
00735 }
00736 }
00737 }
00738
00739 return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
00740 }
00741
00742
00743
00744 }
00745
00746
00747
00748 #endif // OPM_AUTODIFFBLOCK_HEADER_INCLUDED