avoid multiplication with empty matrices.

This commit is contained in:
Robert K 2014-12-03 15:45:49 +01:00
parent 3f821f1d5f
commit 21a9a7c446

View File

@ -30,6 +30,7 @@
#include <vector> #include <vector>
#include <cassert> #include <cassert>
#include <iostream>
namespace Opm namespace Opm
{ {
@ -293,8 +294,18 @@ namespace Opm
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
jac[block] = M( D2.rows(), jac_[block].cols() );
}
else if( jac_[block].nonZeros() == 0 )
jac[block] = D1*rhs.jac_[block];
else if ( rhs.jac_[block].nonZeros() == 0 ) {
jac[block] = D2*jac_[block];
}
else {
jac[block] = D2*jac_[block] + D1*rhs.jac_[block]; jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
} }
}
return function(val_ * rhs.val_, jac); return function(val_ * rhs.val_, jac);
} }
@ -320,8 +331,21 @@ namespace Opm
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
if( jac_[block].nonZeros() == 0 && rhs.jac_[block].nonZeros() == 0 ) {
jac[block] = M( D3.rows(), jac_[block].cols() );
}
else if( jac_[block].nonZeros() == 0 ) {
jac[block] = D3 * ( D1*rhs.jac_[block]);
jac[block] *= -1.0;
}
else if ( rhs.jac_[block].nonZeros() == 0 )
{
jac[block] = D3 * (D2*jac_[block]);
}
else {
jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]); jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
} }
}
return function(val_ / rhs.val_, jac); return function(val_ / rhs.val_, jac);
} }