Alessio Quaglino wrote:
Anyway, I just have a small (but important) bug with uBlas and dolfin
0.6.3 that Johan Jansson couldnt solve so if you (or someone else if you
don't have time) could help me I would really appreciate.
Is the bug present in the later release (0.6.4)?
Garth
Yes. Gbd gives the following output:
Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread -1209726752 (LWP 3498)]
0xb7cd53dc in
boost::numeric::ublas::matrix_assign<boost::numeric::ublas::scalar_assign,
boost::numeric::ublas::basic_full<unsigned int>,
boost::numeric::ublas::generalized_vector_of_vector<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>,
boost::numeric::ublas::vector<boost::numeric::ublas::compressed_vector<double,
0u, boost::numeric::ublas::unbounded_array<unsigned int,
std::allocator<unsigned int> >,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >,
boost::numeric::ublas::unbounded_array<boost::numeric::ublas::compressed_vector<double,
0u, boost::numeric::ublas::unbounded_array<unsigned int,
std::allocator<unsigned int> >,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >,
std::allocator<boost::numeric::ublas::compressed_vector<double, 0u,
boost::numeric::ublas::unbounded_array<unsigned int,
std::allocator<unsigned int> >,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >
, boost::numeric::ublas::compressed_matrix<double,
boost::numeric::ublas::basic_row_major<unsigned int, int>, 0u,
boost::numeric::ublas::unbounded_array<unsigned int,
std::allocator<unsigned int> >,
boost::numeric::ublas::unbounded_array<double, std::allocator<double> > >
(m=@0x8242b7c, e=@0x8242b44) at vector_sparse.hpp:957
957 *it = k_based (i);
The simple example uses the Python interface and the Poisson equation.
I've inserted inside uBlasMatrix.h the following function, which does
exactly what add() does, but with different arguments.
template <class Mat>
inline void uBlasMatrix<Mat>::diagSum(const uBlasVector& x)
{
if( assembled )
{
Assembly_matrix.assign(*this);
assembled = false;
}
else
dolfin_error("Matrix has not been assembled. Did you forget to
call
A.apply()?");
uint s = Assembly_matrix.size1();
uint r = Assembly_matrix.size2();
if( s != r )
dolfin_error("Matrix is not square!");
for( uint i = 0; i < s; i++ )
{
Assembly_matrix(i , i) += x[i];
}
}
My simple test case just import the form file already defined for Poisson
euqation, uses dolfin mesh generator and then assemble the matrix A:
# Create a mesh of the unit square
mesh = UnitSquare(100, 100)
# Import forms (compiled just-in-time with FFC)
forms = import_formfile("Poisson.form")
import poisson as forms
# Assemble linear system
a = forms.PoissonBilinearForm()
A = Matrix()
FEM_assemble(a, A, mesh)
# Sum a vector
b = Vector(A.size(0))
b.operator=(2.0)
A.diagSum(b)
I've also already tried to call A.apply() before but it does not help. The
segmentation fault comes from the following line (inside uBlasMatrix.h):
Assembly_matrix.assign(*this);
However, there is no segmentation fault if I call A.zero() after
FEM_assemble(), so my guess is that this depends on the matrix
representation. Thank you for the help.
Regards,
Alessio Quaglino
Best regards,
Alessio Quaglino
