FullMatrix to SparseMatrix

80 views
Skip to first unread message

Nikki Holtzer

unread,
Oct 7, 2020, 4:51:33 PM10/7/20
to deal.II User Group

Hello everyone,

I am trying to recast a current FullMatrix<double> as a SparseMatrix<double> in order to add this matrix to another SparseMatrix<double>.

I have tried something like this:

SparseMatrix<double> Mat_Sparse;

Mat_Sparse.copy_from(Mat_Full);

Mat_Sparse_2.add(1.,Mat_Sparse);


When doing so I receive the following error:


An error occurred in line <195> of file <dealii-9.1.1/include/deal.II/lac/sparse_matrix.templates.h> in function

    dealii::SparseMatrix<Number>& dealii::SparseMatrix<number>::operator=(double) [with number = double]

The violated condition was:

    cols != nullptr

Additional information:

    (none)



How should I amend this recasting?

Thank you!


Wolfgang Bangerth

unread,
Oct 7, 2020, 4:58:45 PM10/7/20
to dea...@googlegroups.com
On 10/7/20 2:51 PM, Nikki Holtzer wrote:
>
> When doing so I receive the following error:
>
>
> An error occurred in line <195> of file
> <dealii-9.1.1/include/deal.II/lac/sparse_matrix.templates.h> in function
>
>     dealii::SparseMatrix<Number>&
> dealii::SparseMatrix<number>::operator=(double) [with number = double]
>
> The violated condition was:
>
>     cols != nullptr
>
> Additional information:
>
>     (none)

You haven't associated a sparsity pattern with the sparse matrix yet. You'll
have to set that first, and you need to make sure that it has an entry for
each nonzero in the full matrix you are hoping to copy over!

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

Nikki Holtzer

unread,
Oct 8, 2020, 9:16:24 AM10/8/20
to deal.II User Group
When I do the above, with the modifications you have provided, I still obtain the following:

error: no match for ‘operator+’ (operand types are ‘const dealii::SparseMatrix<double>’ and ‘dealii::SparseMatrix<double>’)

or

error: no matching function for call to ‘dealii::SparseMatrix<double>::add(double, dealii::SparseMatrix<double>&) const’


I have not added the 'const' qualifier to my definitions of either of my sparse matrices. 


Nikki Holtzer

unread,
Oct 8, 2020, 5:16:08 PM10/8/20
to deal.II User Group
Hello everyone,

I was just curious if anyone had any suggestions for the above 'const' qualifier problem I have stated above?

Thanks!
Nikki 

Wolfgang Bangerth

unread,
Oct 8, 2020, 9:22:38 PM10/8/20
to dea...@googlegroups.com
On 10/8/20 7:16 AM, Nikki Holtzer wrote:
>
> error: no match for ‘operator+’ (operand types are ‘const
> dealii::SparseMatrix<double>’ and ‘dealii::SparseMatrix<double>’)

SparseMatrix has no operator+, so this will not work.


> error: no matching function for call to
> ‘dealii::SparseMatrix<double>::add(double, dealii::SparseMatrix<double>&) const’

The usual case where this happens is if you are in a member function that is
marked 'const' and the object to which you try to 'add' something is a member
variable. For example in

class X {
SparseMatrix<double> A,B;

void my_function() const // <--- note the 'const'
{
A.add (1.0, B);
}
};

The 'const' on the member function is a promise that you will not change the
values of the members of this class, but A.add(...) will change the member
variable A. The compiler is telling you that this is not allowed.

Nikki Holtzer

unread,
Oct 10, 2020, 5:27:10 PM10/10/20
to deal.II User Group
My sparsity pattern was created in the following way:

DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());

  DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, affine_constraints, true);


which I have used throughout my code with no problems. 


The FullMatrix, Mat1, that I form comes from the outer product of vec1 and vec2 both of size dof_handler.n_dofs. 

I create a Sparse Matrix, Mat2 that i initialize with the sparsity pattern written above. I am not seeing why Mat1 of size dof_handler.n_dofs X dof_handler.n_dofs does not fit into a sparsity pattern that's defined as DynamicSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs());

I am indeed getting an error that there are not enough entries in the SparseMatrix, Mat2.

Thank you!


On Wednesday, October 7, 2020 at 4:58:45 PM UTC-4 Wolfgang Bangerth wrote:

Wolfgang Bangerth

unread,
Oct 11, 2020, 5:38:12 PM10/11/20
to dea...@googlegroups.com

Nikki,
we can continue to speculate, but your description of what is happening is not
specific enough to really know. You'll have to show us the code you use to be
sure.

I don't understand the connection between the code you show and the outer
product. The outer product is a dense matrix, but the sparsity pattern you
create is sparse. Were you hoping to put the outer product of two vectors into
a sparse matrix?

Best
W.
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb8d612f3bf054629da0108d86d6341e7%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637379620364329787&sdata=%2FrnFtYKQyv6SAq3VPanNNWuc44wda22mCApTE3eKs38%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb8d612f3bf054629da0108d86d6341e7%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637379620364334764&sdata=nj5OkstMbznZWPpEdCGgBDRHF1h1tzib22BmHu4Cj7s%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/56cca082-6ae2-4b51-af0d-85ed95661ac4n%40googlegroups.com
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F56cca082-6ae2-4b51-af0d-85ed95661ac4n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb8d612f3bf054629da0108d86d6341e7%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637379620364334764&sdata=xb9XCcWPhiW116ZH7j13ZnGpkwxxCDLVDnLI8pxoBzk%3D&reserved=0>.

Nikki Holtzer

unread,
Oct 11, 2020, 8:09:27 PM10/11/20
to deal.II User Group

I have copied some code below. 


I was hoping to take the outer product which is formed in a FullMatrix, called Kronecker_matrix below, multiply it by a coefficient, and add a matrix to it called nonlinear_part_full. This entire operation I have tried to put into a Sparse matrix, called Nonlinear_part, because I need to add this term to linear_part which is Sparse. I have done this simply because it was the only way I could make the addition between linear_part (Sparse) and nonlinear_part_full (FullMatrix) work all while trying to make the code as efficient as possible. 



template <int dim>

void TimeStep<dim>::prepare()

{ 

  std::cout << "TimeStep<dim>::prepare()" << std::endl;

  const auto &offline_data = *p_offline_data;

  linear_part.reinit(offline_data.sparsity_pattern);

  auto &M_c = offline_data.mass_matrix;

  auto &S_c = offline_data.stiffness_matrix;

  

  linear_part.copy_from(M_c);

  linear_part.add((1. - theta) * kappa, S_c);

  linear_part_inverse.initialize(linear_part);

  nonlinear_part.reinit(offline_data.sparsity_pattern);

  auto &X = offline_data.inner_product_matrix;

  nonlinear_part.copy_from(X);

}


template <int dim>

void TimeStep<dim>::step(Vector<double> &old_solution,double new_t)


{

  const auto &offline_data = *p_offline_data;

  const auto &coefficients = *offline_data.p_coefficients;

  const auto &affine_constraints = offline_data.affine_constraints;

  auto &manufactured = *p_manufactured;


  const auto lam = coefficients.lam;

  std::cout<<"lambda in Time step  = "<<lam<<std::endl;


  const auto M_c = linear_operator(offline_data.mass_matrix);

  const auto S_c = linear_operator(offline_data.stiffness_matrix);

  const auto M_u = linear_operator(offline_data.mass_matrix_unconstrained);

  const auto S_u = linear_operator(offline_data.stiffness_matrix_unconstrained);

  const auto X   = linear_operator(offline_data.inner_product_matrix);


  GrowingVectorMemory<Vector<double>> vector_memory;

  typename VectorMemory<Vector<double>>::Pointer p_new_solution(vector_memory);

  auto &new_solution = *p_new_solution;

  new_solution = old_solution;

  apply_boundary_values(offline_data, coefficients, new_t, new_solution, old_solution);

  

  for (unsigned int m = 0; m < nonlinear_solver_limit; ++m) { Vector<double> residual = M_u * (new_solution - old_solution) 

          + kappa * (1. - theta) * S_u * new_solution +  theta * kappa * S_u * old_solution - kappa * (1.- theta) * lam * new_solution.norm_sqr()                    *  new_solution - kappa * theta * lam * old_solution.norm_sqr() * old_solution;


    affine_constraints.set_zero(residual);

    if (residual.linfty_norm() < nonlinear_solver_tol)

      break;

    double old_solution_normsq = -old_solution.norm_sqr();

     const unsigned int k = old_solution.size();

    FullMatrix<double>Kronecker_matrix(k,k);

    for(unsigned int i=0; i<k; ++i)

        for(unsigned int j=0; j<k; ++j)

                Kronecker_matrix(i,j) = old_solution[i]*old_solution[j];


    FullMatrix<double>nonlinear_part_full;

    nonlinear_part_full.copy_from(nonlinear_part);

    double nonlinear_part_coef = - kappa * (1.- theta) * lam * old_solution_normsq;

    double Kronecker_matrix_coef = -2.* lam * kappa * (1- theta);

    nonlinear_part_full *= nonlinear_part_coef; // -k(1- theta) lambda |Psi|^2*X

   Kronecker_matrix *= Kronecker_matrix_coef; //(-2.* lam * kappa * (1- theta)* Kronecker_matrix)

   nonlinear_part_full.add(1.,Kronecker_matrix); //-k(1- theta) lambda |Psi|^2*X - 2k(1-theta)lambda Kronecker_matrix

   

   SparseMatrix<double> Nonlinear_part;

    Nonlinear_part.reinit(offline_data.sparsity_pattern);

    Nonlinear_part.copy_from(nonlinear_part_full);

    linear_part.add(1., Nonlinear_part);

    auto system_matrix = linear_operator(linear_part);



    SolverControl solver_control(linear_solver_limit, linear_solver_tol);

    SolverGMRES<> solver(solver_control);


    auto system_matrix_inverse = inverse_operator(system_matrix, solver, linear_part_inverse);

    Vector<double> update = system_matrix_inverse * (-1. * residual);


    affine_constraints.set_zero(update);



    new_solution += update;

Wolfgang Bangerth

unread,
Oct 12, 2020, 7:27:35 PM10/12/20
to dea...@googlegroups.com
On 10/11/20 6:09 PM, Nikki Holtzer wrote:
> I was hoping to take the outer product which is formed in a FullMatrix, called
> Kronecker_matrix below, multiply it by a coefficient, and add a matrix to it
> called nonlinear_part_full. This entire operation I have tried to put into a
> Sparse matrix, called Nonlinear_part, because I need to add this term to
> linear_part which is Sparse. I have done this simply because it was the only
> way I could make the addition between linear_part (Sparse) and
> nonlinear_part_full (FullMatrix) work all while trying to make the code as
> efficient as possible.

But if the nonlinear part is a dense matrix, you will necessarily add entries
into the sparse matrix that are not present in the sparsity pattern. When you
build a sparse matrix with a certain sparsity pattern, you are saying "these
are the entries that could possibly be nonzero, please allocate memory for
them". So when you later add the dense matrix to the sparse one, you're
breaking this promise.

I don't quite understand why you copy everything into a sparse matrix in your
code. Couldn't you just keep everything in dense matrices and do the linear
solve with that?

Nikki Holtzer

unread,
Oct 13, 2020, 12:28:07 PM10/13/20
to deal.II User Group
The issue that I have when  all matrices are FullMatrix matrices is I can't seem to get the .add function to result in anything other than error. If I keep ONLY both nonlinear_part and linear_part as FullMatrix 
I receive:

error: no matching function for call to ‘dealii::FullMatrix<double>::add(double, const dealii::SparseMatrix<double>&)’

   linear_part.add((1. - theta) * kappa, S_c);


Which leads me to believe that I cannot add a SparseMatrix (S_c) to a FullMatrix (linear_part). I tried to back track and change S_c to be a FullMatrix instead of a SparseMatrix 

For instance: 


FullMatrix<double> mass_matrix_unconstrainted;
mass_matrix_unconstrained = 0;
FullMatrix<double> cell_mass_matrix(dofs_per_cell, dofs_per_cell);
cell_mass_matrix = 0;


Then I fill cell_mass_matrix.

mass_matrix_unconstrained.add(local_dof_indices, cell_mass_matrix);

error: no matching function for call to ‘dealii::FullMatrix<double>::add(std::vector<unsigned int>&, dealii::FullMatrix<double>&)’

     mass_matrix_unconstrained.add(local_dof_indices, cell_mass_matrix);


Alternatively if I try:

mass_matrix_unconstrained = mass_matrix_unconstrained  + local_dof_indices*cell_mass_matrix;


I receive.

error: no match for ‘operator*’ (operand types are ‘std::vector<unsigned int>’ and ‘dealii::FullMatrix<double>’)

     mass_matrix_unconstrained = mass_matrix_unconstrained + local_dof_indices * cell_mass_matrix;


I thought maybe this is a type error, i.e. I can't multiply a vector of int by a matrix of double but I can't seem to be able to convert the vector without error either. 


Thank you,

Nikki



Wolfgang Bangerth

unread,
Oct 13, 2020, 3:37:03 PM10/13/20
to dea...@googlegroups.com
On 10/13/20 10:28 AM, Nikki Holtzer wrote:
> The issue that I have when  all matrices are FullMatrix matrices is I can't
> seem to get the .add function to result in anything other than error. If I
> keep ONLY both nonlinear_part and linear_part as FullMatrix
> I receive:
>
> error: no matching function for call to
> ‘dealii::FullMatrix<double>::add(double, const dealii::SparseMatrix<double>&)’
>
>    linear_part.add((1. - theta) * kappa, S_c);
>
>
> Which leads me to believe that I cannot add a SparseMatrix (S_c) to a
> FullMatrix (linear_part).

Correct -- there is no such function.


> I tried to back track and change S_c to be a
> FullMatrix instead of a SparseMatrix
>
> For instance:
>
>
> FullMatrix<double> mass_matrix_unconstrainted;
> mass_matrix_unconstrained = 0;
> FullMatrix<double> cell_mass_matrix(dofs_per_cell, dofs_per_cell);
> cell_mass_matrix = 0;
>
>
> Then I fill cell_mass_matrix.
>
> mass_matrix_unconstrained.add(local_dof_indices, cell_mass_matrix);
>
> error: no matching function for call to
> ‘dealii::FullMatrix<double>::add(std::vector<unsigned int>&,
> dealii::FullMatrix<double>&)’
>
>      mass_matrix_unconstrained.add(local_dof_indices, cell_mass_matrix);

Right, that function also doesn't exist, but you could do
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
mass_matrix_unconstrainted(local_dof_indices[i], local_dof_indices[j])
+= cell_mass_matrix(i,j);


> Alternatively if I try:
>
> mass_matrix_unconstrained = mass_matrix_unconstrained  +
> local_dof_indices*cell_mass_matrix;
>
>
> I receive.
>
> error: no match for ‘operator*’ (operand types are ‘std::vector<unsigned int>’
> and ‘dealii::FullMatrix<double>’)

Right. There is also no operator* for vector<unsigned int> and
FullMatrix<double>. It's also unclear to me what exactly that would actually
supposed to do.

Writing the double-loop above is the way to go.

Nikki Holtzer

unread,
Oct 14, 2020, 10:14:34 AM10/14/20
to deal.II User Group

Ok thank you that has resolved most of the adding issue. I am still receiving the following error. It reads to me that there's still an error on the 'FullMatrix::add()' function. I have verified that they are all gone however. What else should I be looking for?




An error occurred in line <1335> of file </home/u7/nholtzer/dealii-9.1.1/include/deal.II/lac/full_matrix.h> in function

    void dealii::FullMatrix< <template-parameter-1-1> >::add(dealii::FullMatrix< <template-parameter-1-1> >::size_type, dealii::FullMatrix< <template-parameter-1-1> >::size_type, const index_type*, const number2*, bool, bool) [with number2 = double; index_type = unsigned int; number = double; dealii::FullMatrix< <template-parameter-1-1> >::size_type = long unsigned int]

The violated condition was:

    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(row), decltype(this->m())>::type)>:: type>(row) < static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(row), decltype(this->m())>::type)>:: type>(this->m())

Additional information:

    Index 0 is not in the half-open range [0,0). In the current case, this half-open range is in fact empty, suggesting that you are accessing an element of an empty collection such as a vector that has not been set to the correct size.


Stacktrace:

-----------

#0  /home/u7/nholtzer/dealii-9.1.1/lib/libdeal_II.g.so.9.1.1: void dealii::FullMatrix<double>::add<double, unsigned int>(unsigned long, unsigned long, unsigned int const*, double const*, bool, bool)

#1  /home/u7/nholtzer/dealii-9.1.1/lib/libdeal_II.g.so.9.1.1: void dealii::AffineConstraints<double>::distribute_local_to_global<dealii::FullMatrix<double>, dealii::Vector<double> >(dealii::FullMatrix<double> const&, dealii::Vector<double> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::FullMatrix<double>&, dealii::Vector<double>&, bool, std::integral_constant<bool, false>) const

#2  ./gravwave: void dealii::AffineConstraints<double>::distribute_local_to_global<dealii::FullMatrix<double> >(dealii::FullMatrix<double> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::FullMatrix<double>&) const

#3  ./gravwave: OfflineData<1>::assemble()

#4  ./gravwave: OfflineData<1>::prepare()

#5  ./gravwave: main

Wolfgang Bangerth

unread,
Oct 14, 2020, 10:31:51 AM10/14/20
to dea...@googlegroups.com
On 10/14/20 8:14 AM, Nikki Holtzer wrote:
>     void dealii::FullMatrix< <template-parameter-1-1>
> >::add(dealii::FullMatrix< <template-parameter-1-1> >::size_type,
> dealii::FullMatrix< <template-parameter-1-1> >::size_type, const index_type*,
> const number2*, bool, bool) [with number2 = double; index_type = unsigned int;
> number = double; dealii::FullMatrix< <template-parameter-1-1> >::size_type =
> long unsigned int]
>
> The violated condition was:
>
>     static_cast<typename ::dealii::internal::argument_type<void( typename
> std::common_type<decltype(row), decltype(this->m())>::type)>:: type>(row) <
> static_cast<typename ::dealii::internal::argument_type<void( typename
> std::common_type<decltype(row), decltype(this->m())>::type)>:: type>(this->m())
>
> Additional information:
>
>     Index 0 is not in the half-open range [0,0). In the current case, this
> half-open range is in fact empty, suggesting that you are accessing an element
> of an empty collection such as a vector that has not been set to the correct size.

Reading the error message, it seems to me that you didn't size the matrix
correctly.

Nikki Holtzer

unread,
Oct 14, 2020, 8:14:46 PM10/14/20
to deal.II User Group
The first error seems to be occurring at:

affine_constraints.distribute_local_to_global(cell_nodal_mass_matrix, local_dof_indices, nodal_mass_matrix);

I have verified that cell_nodal_mass_matrix is indeed the same size as nodal_mass_matrix.

cell_nodal_mass_matrix is initialized as:


FullMatrix<double> cell_nodal_mass_matrix(dofs_per_cell, dofs_per_cell);


and nodal_mass_matrix has the same dimensions. 

local_dof_indices seem to be the problem. I see that on dealii 's website that  as long as cell_nodal_mass_matrix and local_dof_indices have the same number of elements this should work. 

I have defined local_dof_indices in the following way:


std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

fe_values.reinit(cell);

cell->get_dof_indices(local_dof_indices);


The number of elements are:

cell_nodal_mass_matrix: 136

local_dof_indices: 24

How would I resolve this?

Thank you,

Nikki 

Wolfgang Bangerth

unread,
Oct 16, 2020, 10:37:40 AM10/16/20
to dea...@googlegroups.com

Nikki,
I can't tell from these pieces of code what dimension is wrong. But you will
find it very useful to run your program in a debugger so that you can see
which matrix access it is that triggers the error, and why. There are a number
of video lectures about running a program in a debugger.

Best
WB
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0f4b5d9bc81e4aa9984708d8709f5570%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637383176932210123&sdata=C2kpiS1UGPD2KcgJ6OUXv08tykRiVvHFCtyfcGgWXmE%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0f4b5d9bc81e4aa9984708d8709f5570%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637383176932210123&sdata=07ZJ%2FIjSK0cZXfjqTRoqRNN6HVktn%2FPCDgvP20HbdnQ%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/a9c2f25b-e5fc-4fe1-92e7-9641331a6143n%40googlegroups.com
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fa9c2f25b-e5fc-4fe1-92e7-9641331a6143n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0f4b5d9bc81e4aa9984708d8709f5570%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637383176932220121&sdata=ICBH4g%2FL1QvA0hylXhcv6QKWzZ5OTg1YgQZckIEOHxU%3D&reserved=0>.

Nikki Holtzer

unread,
Oct 16, 2020, 12:35:25 PM10/16/20
to deal.II User Group
I was able to diagnose and fix the problem. Thank you. I still have one remaining error since changing all my matrices to FullMatrix. I have diagnosed that the issue is the very last line of code below, when attempting to update my system. I have printed the line in red.

for (unsigned int m = 0; m < nonlinear_solver_limit; ++m) {

    Vector<double> residual = M_u * (new_solution - old_solution) + kappa * (1. - theta) * S_u * new_solution                                                       +   theta * kappa * S_u * old_solution - kappa * (1.- theta) * lam * new_solution.norm_sqr() * new_solution - kappa * theta * lam * old_solution.norm_sqr() * old_solution;

    affine_constraints.set_zero(residual);


    if (residual.linfty_norm() < nonlinear_solver_tol)

      break;

    double new_solution_normsq = -new_solution.norm_sqr();

    const unsigned int k = new_solution.size();

    FullMatrix<double>Kronecker_matrix(k,k);

    for(unsigned int i=0; i<k; ++i)

        for(unsigned int j=0; j<k; ++j)

        Kronecker_matrix(i,j) = new_solution[i]*new_solution[j];

    FullMatrix<double>nonlinear_part_full;

    nonlinear_part_full.copy_from(nonlinear_part);

    FullMatrix<double>nonlinear_part_inverse;

    double nonlinear_part_coef = - kappa * (1.- theta) * lam * new_solution_normsq;

    double Kronecker_matrix_coef = -2.* lam * kappa * (1- theta);

    nonlinear_part_full *= nonlinear_part_coef; // -k(1- theta) lambda |Psi|^2*X

    Kronecker_matrix *= Kronecker_matrix_coef; //(-2.* lam * kappa * (1- theta)* Kronecker_matrix)

    for(unsigned int i=0; i<dofs_per_cell; ++i)

        for(unsigned int j=0; j<dofs_per_cell; ++j)

                nonlinear_part_full(i,j) += Kronecker_matrix(i,j);

    for(unsigned int i=0; i<dofs_per_cell; ++i)

        for(unsigned int j=0; j<dofs_per_cell; ++j)

                nonlinear_part_full(i,j) += linear_part(i,j);

    auto system_matrix = linear_operator(nonlinear_part_full);

   SolverControl solver_control(linear_solver_limit, linear_solver_tol);

   SolverGMRES<> solver(solver_control);


    auto system_matrix_inverse = inverse_operator(system_matrix, solver, nonlinear_part_inverse);

    Vector<double> update = system_matrix_inverse * (-1. * residual);


The error I'm receiving is:

An error occurred in line <179> of file </home/u7/nholtzer/dealii-9.1.1/include/deal.II/lac/full_matrix.templates.h> in function

 void dealii::FullMatrix<number>::vmult(dealii::Vector<OtherNumber>&, const dealii::Vector<OtherNumber>&, bool) const [with number2 = double; number = double]

The violated condition was:

    !this->empty()

Additional information:

    (none)


I have also tried this approach 

Vector<double> update;
system_matrix_inverse.vmult(update, -1*residual)

but it resulted in the same error. 

Wolfgang Bangerth

unread,
Oct 17, 2020, 2:32:02 PM10/17/20
to dea...@googlegroups.com
On 10/16/20 10:35 AM, Nikki Holtzer wrote:
>
> The error I'm receiving is:
>
> An error occurred in line <179> of file
> </home/u7/nholtzer/dealii-9.1.1/include/deal.II/lac/full_matrix.templates.h>
> in function
>
>  void dealii::FullMatrix<number>::vmult(dealii::Vector<OtherNumber>&, const
> dealii::Vector<OtherNumber>&, bool) const [with number2 = double; number = double]
>
> The violated condition was:
>
>     !this->empty()
>
> Additional information:
>
>     (none)

The error message states that you are trying to multiply a vector by a matrix
that is empty (that is, has zero columns and/or rows). Have you run the code
in a debugger to see which matrix this is?

Nikki Holtzer

unread,
Oct 17, 2020, 4:12:44 PM10/17/20
to deal.II User Group
Yes it appears to be the system_matrix_inverse because I can print out the residual and it is indeed not empty. 

Wolfgang Bangerth

unread,
Oct 18, 2020, 10:40:42 PM10/18/20
to dea...@googlegroups.com
On 10/17/20 2:12 PM, Nikki Holtzer wrote:
> Yes it appears to be the system_matrix_inverse because I can print out the
> residual and it is indeed not empty.

There is little any of can help you with without having access to the actual
code. The usual way to approach this now is to make the program smaller and
smaller (e.g., removing the entire code that actually puts anything into the
matrices) until you end up with something that has as few lines as possible
while still showing the error. You should typically be able to get a program
to ~100 lines. Hopefully, it is easier to see with such a minimal test case
what is wrong.
Reply all
Reply to author
Forward
0 new messages