[deal.II] Re: LinearOperator of BlockVector issue and potential fix
Dear all,
The actual fix to this issue is much simpler than what I originally
thought. I just need to make sure the destination BlockVector *v* has the
proper block sizes when I define the vmult(),
vmult_add(),
Tvmult(),
Tvmult_add(),
reinit_range_vector,
and reinit_domain_vector for the LinearOperator.
We don't need to touch the definition of operator*() at all. I attach the
following test case as an example, in which I use two FullMatrix objects to
construct a LinearOperator that operates on BlockVector.
using namespace dealii;
// LinearOperator
using Payload =
dealii::internal::LinearOperatorImplementation::EmptyPayload;
LinearOperator, BlockVector, Payload>
usr_linearOperator_1;
// For demonstration purpose, define LinearOperator using two FullMatrix
FullMatrix usr_matrix_1(3, 3);
usr_matrix_1(0,0) = 1;
usr_matrix_1(0,1) = 2;
usr_matrix_1(0,2) = 3;
usr_matrix_1(1,0) = 4;
usr_matrix_1(1,1) = 5;
usr_matrix_1(1,2) = 6;
usr_matrix_1(2,0) = 7;
usr_matrix_1(2,1) = 8;
usr_matrix_1(2,2) = 9;
FullMatrix usr_matrix_2(2, 2);
usr_matrix_2(0,0) = -1;
usr_matrix_2(0,1) = -2;
usr_matrix_2(1,0) = -3;
usr_matrix_2(1,1) = -4;
usr_linearOperator_1.vmult = [& usr_matrix_1, &
usr_matrix_2](BlockVector , const BlockVector ) {
std::vector sizes_per_block(u.n_blocks());
for (unsigned int i = 0; i < u.n_blocks(); ++i)
sizes_per_block[i] = u.block(i).size();
v.reinit(sizes_per_block);
usr_matrix_1.vmult(v.block(0), u.block(0));
usr_matrix_2.vmult(v.block(1), u.block(1));
};
usr_linearOperator_1.vmult_add =
[_linearOperator_1](BlockVector , const BlockVector
) {
std::vector sizes_per_block(u.n_blocks());
for (unsigned int i = 0; i < u.n_blocks(); ++i)
sizes_per_block[i] = u.block(i).size();
BlockVector temp_v(sizes_per_block);
usr_linearOperator_1.vmult(temp_v, u);
v += temp_v;
};
usr_linearOperator_1.Tvmult = [& usr_matrix_1, &
usr_matrix_2](BlockVector , const BlockVector ) {
std::vector sizes_per_block(u.n_blocks());
for (unsigned int i = 0; i < u.n_blocks(); ++i)
sizes_per_block[i] = u.block(i).size();
v.reinit(sizes_per_block);
usr_matrix_1.Tvmult(v.block(0), u.block(0));
usr_matrix_2.Tvmult(v.block(1), u.block(1));
};
usr_linearOperator_1.Tvmult_add =
[usr_linearOperator_1](BlockVector , const BlockVector
) {
std::vector sizes_per_block(u.n_blocks());
for (unsigned int i = 0; i < u.n_blocks(); ++i)
sizes_per_block[i] = u.block(i).size();
BlockVector temp_v(sizes_per_block);
usr_linearOperator_1.Tvmult(temp_v, u);
v += temp_v;
};
usr_linearOperator_1.reinit_range_vector = [& usr_matrix_1, &
usr_matrix_2](BlockVector , bool omit_zeroing_entries) {
v.reinit(2);
(v.block(0)).reinit(usr_matrix_1.m(), omit_zeroing_entries);
(v.block(1)).reinit(usr_matrix_2.m(), omit_zeroing_entries);
};
usr_linearOperator_1.reinit_domain_vector = [& usr_matrix_1, &
usr_matrix_2](BlockVector , bool omit_zeroing_entries) {
v.reinit(2);
(v.block(0)).reinit(usr_matrix_1.n(), omit_zeroing_entries);
(v.block(1)).reinit(usr_matrix_2.n(), omit_zeroing_entries);
};
Best,
Tao
On Saturday, April 27, 2024 at 10:33:11 AM UTC-4 Tao Jin wrote:
> Dear all,
>
> *Background*: I have a BlockVector *u* that needs to be updated
> via a series of rank-2 update to get BlockVector *v*. For
> notation convenience, let's call this series of rank-2 update as *P*. I
> can define *P *as either a LinearOperator or BlockLinearOperator, since I
> need to solve a linear system later that involves *P* using an iterative
> solver (for example, conjugate-gradient method). Defining *P *as a
> BlockLinearOperator does not work in this case, since I only know how to
> perform the rank-2 updates on *u* but not the block structure of the
> operator *P*.
>
> *Solution*: My solution is to define *P* as a LinearOperator that
> operates on the source BlockVector *u* to get the destination BlockVector
> *v*. I know this is unconventional, because LinearOperator usually
> operates on Vector but not BlockVector. But theoretically, a LinearOperator
> should still be able to operator on a BlockVector to produce another
> BlockVector.
>
> *Issue*: I can define the needed LinearOperator *P1 *and* P2 *that
> operate on a BlockVector without a problem. I can call *P1*.vmult(*v*, *u*)
> and *P2*.vmult(*v*, *u*) correctly. However, when I define a composite
> operator such as (for demonstration purpose)
> *auto total_op = P1 * P2; *
> and then call
> *total_op.vmult(v, u);*
> I get the following error message:
> """
>
> An error occurred in line <1487> of file
> in
> function
> dealii::BlockVectorBase::BlockType&
[deal.II] LinearOperator of BlockVector issue and potential fix
Dear all,
*Background*: I have a BlockVector *u* that needs to be updated via
a series of rank-2 update to get BlockVector *v*. For notation
convenience, let's call this series of rank-2 update as *P*. I can define *P
*as either a LinearOperator or BlockLinearOperator, since I need to solve a
linear system later that involves *P* using an iterative solver (for
example, conjugate-gradient method). Defining *P *as a BlockLinearOperator
does not work in this case, since I only know how to perform the rank-2
updates on *u* but not the block structure of the operator *P*.
*Solution*: My solution is to define *P* as a LinearOperator that operates
on the source BlockVector *u* to get the destination BlockVector *v*. I
know this is unconventional, because LinearOperator usually operates on
Vector but not BlockVector. But theoretically, a LinearOperator should
still be able to operator on a BlockVector to produce another BlockVector.
*Issue*: I can define the needed LinearOperator *P1 *and* P2 *that operate
on a BlockVector without a problem. I can call *P1*.vmult(*v*, *u*) and *P2*
.vmult(*v*, *u*) correctly. However, when I define a composite operator
such as (for demonstration purpose)
*auto total_op = P1 * P2; *
and then call
*total_op.vmult(v, u);*
I get the following error message:
"""
An error occurred in line <1487> of file
in
function
dealii::BlockVectorBase::BlockType&
dealii::BlockVectorBase::block(unsigned int) [with VectorType =
dealii::Vector; dealii::BlockVectorBase::BlockType =
dealii::Vector]
The violated condition was:
::dealii::deal_II_exceptions::internals::compare_less_than(i,
n_blocks())
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.
"""
*Fix*: This error message is due to the definition of operator* (Line 582
in linear_operator.h)
template
LinearOperator
operator*(const LinearOperator & first_op,
const LinearOperator _op)
{
..
return_op.vmult = [first_op, second_op](Range , const Domain ) {
GrowingVectorMemory vector_memory;
typename VectorMemory<*Intermediate*>::Pointer i(vector_memory);
second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries
=*/true);
second_op.vmult(*i, u);
first_op.vmult(v, *i);
..
}
It looks like that *Intermediate* only has the correct memory size but not
the block information. My fix is to add one line of code
*(*i).reinit(u.n_blocks());*
in my own definition operator*(). This fixes the error and works correctly.
*Questions*:
1. Does it even make sense to define a LinearOperator operating on a
BlockVector? If yes, should we expand the dealii capability to consider
this option?
2. For my rank-2 update operation *P*, I cannot define it as a
BlockLinearOperator. Even though *P *operates on a BlockVector, but *P *itself
does not possess blocks. Any thoughts?
3. Another solution is to convert BlockVector *u* to a Vector and define *P
*as a LinearOperator operating on Vector. After the rank-2 update, convert
the Vector back to a BlockVector. However, it seems rather awkward. Does
the community have better ideas?
I have attached my code and a test case here.
Best,
Tao
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4f0dfa54-cec3-4725-ae17-a2ef155e94a3n%40googlegroups.com.
/* -
*
* Copyright (C) 2019 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* -
* This tutorial program was contributed by Martin Kronbichler
*/
#include
#include
#include
#include
#include
#include
#include
using namespace dealii;
// define my own operator*() to make the following line work:
// const auto total_op = op_1 * op_2;
// total_op.vmult(v, u);
// where v and u are BlockVector, and op_1 and op_2 are LinearOperator that
Re: [deal.II] LinearOperator for fully dense and non-square matrix
Dear Wolfgang, I am glad to report that the LinearOperator modules also work for fully dense and non-rectangular matrices (FullMatrix). I can now use the inverse operator to solve the linear system: B = B_0 + W * M * W^T. B_0 is an n by n sparse symmetric positive-definite matrix (n is large > 100k); M is a m by m *full matrix* (m is small m = 20 for example); W is a n by m full matrix; Thanks again for your advice. Best, Tao On Saturday, April 20, 2024 at 12:27:02 AM UTC-4 Tao Jin wrote: > Dear Wolfgang, > > Thanks for the advice! > Best, > > Tao > > On Friday, April 19, 2024 at 11:58:33 PM UTC-4 Wolfgang Bangerth wrote: > >> On 4/19/24 21:54, Tao Jin wrote: >> > >> > Now that my LinearOperator is defined as B_0 + W * M * W^T, I only need >> to >> > find an appropriate iterative solver that can solve the inverse of this >> > operator (the underlying matrix is fully dense) efficiently. >> >> It is entirely unimportant that the matrix is dense. For iterative >> solvers, >> this is of no consequence at all. What matters is whether, for example, >> it is >> symmetric and/or positive definite. The other things that matters is >> whether >> you can construct a good preconditioner for it. >> >> You really should read through step-20 and step-22 to see how all of this >> works in practice. >> >> Best >> W. >> >> -- >> >> Wolfgang Bangerth email: bang...@colostate.edu >> www: http://www.math.colostate.edu/~bangerth/ >> >> >> -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/cce59280-c6f2-428b-97ce-9335b78bf90fn%40googlegroups.com.
Re: [deal.II] LinearOperator for fully dense and non-square matrix
Dear Wolfgang, Thanks for the advice! Best, Tao On Friday, April 19, 2024 at 11:58:33 PM UTC-4 Wolfgang Bangerth wrote: > On 4/19/24 21:54, Tao Jin wrote: > > > > Now that my LinearOperator is defined as B_0 + W * M * W^T, I only need > to > > find an appropriate iterative solver that can solve the inverse of this > > operator (the underlying matrix is fully dense) efficiently. > > It is entirely unimportant that the matrix is dense. For iterative > solvers, > this is of no consequence at all. What matters is whether, for example, it > is > symmetric and/or positive definite. The other things that matters is > whether > you can construct a good preconditioner for it. > > You really should read through step-20 and step-22 to see how all of this > works in practice. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/544b149c-aab3-4bc1-a11a-0934ea8d4517n%40googlegroups.com.
Re: [deal.II] LinearOperator for fully dense and non-square matrix
Dear Wolfgang, What you said makes sense. Now that my LinearOperator is defined as B_0 + W * M * W^T, I only need to find an appropriate iterative solver that can solve the inverse of this operator (the underlying matrix is fully dense) efficiently. Best, Tao On Friday, April 19, 2024 at 11:25:10 PM UTC-4 Wolfgang Bangerth wrote: > On 4/19/24 19:42, Tao Jin wrote: > > > > Thank you so much for your reply. I think my problem has two potential > issues: > > 1. W * M * W^T will be an n by n fully dense matrix with n > 100k. Even > if it > > is valid to use LinearOperator to perform W * M * W^T, will memory > required to > > store this operator be an issue? > > 2. I still have to solve B x= (B_0 + W * M * W^T) x = p. Whether linear > > solvers such as CG will be efficient enough is also an issue since W * M > * W^T > > is fully dense. > > Tao, > like I mentioned, a line like > S = linear_operator(W) * linear_operator(M) * transpose_operator(W) > does not actually compute the product of these matrices. The documentation > of > step-20 and step-22 goes to great length in explaining this issue. All the > object S provides is the ability to multiple S by a vector, for which you > only > need matrix-vector products because > S x > = (W M W^T) x > = W (M (W^T x)) > > The product of matrices is never computed because it is never needed. > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/c9d2909f-749f-422d-942b-4cb0ca78b6edn%40googlegroups.com.
Re: [deal.II] LinearOperator for fully dense and non-square matrix
Hi Wolfgang, Thank you so much for your reply. I think my problem has two potential issues: 1. W * M * W^T will be an n by n fully dense matrix with n > 100k. Even if it is valid to use LinearOperator to perform W * M * W^T, will memory required to store this operator be an issue? 2. I still have to solve B x= (B_0 + W * M * W^T) x = p. Whether linear solvers such as CG will be efficient enough is also an issue since W * M * W^T is fully dense. Anyhow, I will do some experiments and let the community know. Best, Tao On Friday, April 19, 2024 at 6:58:18 PM UTC-4 Wolfgang Bangerth wrote: > On 4/18/24 12:06, Tao Jin wrote: > > > > I am developing a quasi Newton solver in deal.ii. During each iteration, > I > > need to solve a linear system B * x = p, where the matrix B has the > following > > form: > > B = B_0 + W * M * W^T. > > B_0 is an n by n sparse matrix (n is large > 100k); > > M is a m by m *full matrix* (m is small m = 20 for example); > > W is a n by m full matrix; > > Since W * M * W^T will be an n by n fully dense matrix, I cannot afford > to > > explicitly store it. > > > > Does the deal.ii community have any suggestion what is the most > efficient way > > to solve this linear system? I am thinking to use LinearOperator to > define and > > operate on W * M * W^T, but I am not sure whether it supports full > matrix and > > how efficient it would be. > > Tao: > I'm not an expert in the LinearOperator framework, but operators such as > yours > appear, for example, in step-20: > > https://dealii.org/developer/doxygen/deal.II/step_20.html#step_20-Implementationoflinearsolversandpreconditioners > There, opS has exactly the form you describe, and it does not actually > compute > the entries of the S matrix, but only represents the *action* as three > matrix-vector products. > > That only leaves the question of whether you can use linear_operator(M) on > full matrices M. I don't know the answer to that, but what happens if you > try? > > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/b37ca7e7-1eac-4c23-9e71-9f8d539b07efn%40googlegroups.com.
[deal.II] LinearOperator for fully dense and non-square matrix
Dear all, I am developing a quasi Newton solver in deal.ii. During each iteration, I need to solve a linear system B * x = p, where the matrix B has the following form: B = B_0 + W * M * W^T. B_0 is an n by n sparse matrix (n is large > 100k); M is a m by m *full matrix* (m is small m = 20 for example); W is a n by m full matrix; Since W * M * W^T will be an n by n fully dense matrix, I cannot afford to explicitly store it. Does the deal.ii community have any suggestion what is the most efficient way to solve this linear system? I am thinking to use LinearOperator to define and operate on W * M * W^T, but I am not sure whether it supports full matrix and how efficient it would be. Best, Tao -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/b7aed966-09b5-4486-ba68-3d74183fc203n%40googlegroups.com.
Re: [deal.II] Solution transfer for BlockVector
Sorry, my bad. It works for BlockVector. However, there is a discrepancy on deal.ii manual. Below is from the deal.ii manual about SolutionTransfer module: https://www.dealii.org/current/doxygen/deal.II/classSolutionTransfer.html#a1171abf820e35e02304017ef3daeeaa5 """ Although the refine_interpolate functions are allowed to be called multiple times, e.g. for interpolating several solution vectors, there is the following possibility of interpolating several functions simultaneously. std::vector > solutions_old(n_vectors, Vector <https://www.dealii.org/current/doxygen/deal.II/classVector.html> (n)); ... std::vector <https://www.dealii.org/current/doxygen/deal.II/classVector.html> > solutions(n_vectors, Vector <https://www.dealii.org/current/doxygen/deal.II/classVector.html> (n)); *soltrans.refine_interpolate(solutions_old, solutions);* """ However, *refine_interpolate() cannot take *std::vector. Only *interpolate*() has an interface for std::vector. Best, Tao On Tuesday, January 23, 2024 at 4:23:15 PM UTC-5 Wolfgang Bangerth wrote: > On 1/23/24 13:03, Tao Jin wrote: > > > > I am trying to use the SolutionTransfer capability of dealii to > interpolate a > > solution from an old mesh to a new mesh. The solution is stored in a > > BlockVector. My understanding is that current solution transfer > functions can > > only interpolate Vector, not BlockVector. Am I right? > > No, this is not true. Have you tried? It should work just fine. > Best > W. > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/df3690c7-1327-429d-864f-2745e1339a90n%40googlegroups.com.
[deal.II] Solution transfer for BlockVector
Dear all, I am trying to use the SolutionTransfer capability of dealii to interpolate a solution from an old mesh to a new mesh. The solution is stored in a BlockVector. My understanding is that current solution transfer functions can only interpolate Vector, not BlockVector. Am I right? Best, Tao -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/e9cef0a5-e192-4ff7-b493-f3dd88a6cd23n%40googlegroups.com.
Re: [deal.II] LAPACK error in syevx
Hi Wolfgang,
Thanks for your input. I would be delighted to look into this issue and
potentially make a contribution.
Also, I noticed that SymmetricTensor<2, dim> also has a function called
*eigenvectors,
*which could robustly provide the eigenvalues and eigenvectors for a
symmetric 2nd order tensor. This function performs well regardless of the
scaling. Below is my function to obtain the eigen system of a 2nd order
symmetric tensor:
template
void spectrum_decomposition(SymmetricTensor<2, dim> const &
symmetric_tensor,
Vector & myEigenvalues,
std::vector> &
myEigenvectors)
{
std::array< std::pair< double, Tensor< 1, dim > >, dim > myEigenSystem;
myEigenSystem = *eigenvectors*(symmetric_tensor);
for (int i = 0; i < dim; i++)
{
myEigenvalues[i] = myEigenSystem[i].first;
myEigenvectors[i] = myEigenSystem[i].second;
}
}
Best,
Tao
On Wednesday, June 14, 2023 at 10:58:37 PM UTC-4 Wolfgang Bangerth wrote:
>
> Tao:
> Yes, this is cumbersome indeed. I looked into this file, and it contains a
> number of places where deal.II outputs to std::cerr. This is not the style
> we
> favor today. I looked it up, and these lines were written in 2005,
> probably
> with the very first version of that file.
>
> I have opened an issue about this:
> https://github.com/dealii/dealii/issues/15363
> As with all things deal.II, we would gladly help any kind of help we can
> get
> -- so if the error message bothers you, perhaps this is your project: Fix
> this
> particular issue, and get it merged into deal.II!
>
> Underlying your question is the issue that dependent on the value of the
> scaling, you either do or do not get an error. I believe that that is
> because
> you need to also scale the upper/lower bounds and the tolerance with which
> you
> call the LAPACK function. When you multiply the 2x2 matrix by a scalar
> factor,
> all of the eigenvalues also scale and consequently the tolerance should
> also
> scale.
>
> Best
> W.
>
> On 6/14/23 09:47, Tao Jin wrote:
> > *** Caution: EXTERNAL Sender ***
> >
> > Dear all,
> >
> > I received the message "LAPACK error in syevx" when using the
> > /*compute_eigenvalues_symmetric*/ of */lapack_full_matrix/* (the
> interface to
> > LAPACK in dealii) to compute the eigenvalues of a symmetric matrix with
> > relatively small matrix entries. I have attached a test file (main.cc)
> in the
> > end.
> >
> > To provide more information, I construct the following 2x2 symmetric
> matrix
> > const double *scalar = 1.0;*
> > for (unsigned int i = 0; i < dim; i++)
> > for (unsigned int j = 0; j <= i; j++)
> > {
> > symmetric_tensor[i][j] = (i + j + 1) * *scalar*;
> > if (j != i)
> > symmetric_tensor[j][i] = symmetric_tensor[i][j];
> > }
> >
> > *In this case, the eigenvalue decomposition runs without error message:*
> > Spectrum decomposition of a symmetric 2nd order tensor
> > Eigenvalues: -0.236068 4.23607
> > Spectrum decomposition succeeded!
> >
> > However, if I set
> > const double *scalar = 1.0e-3; // to make the matrix entries
> relatively small*
> > The eigenvalue decomposition is still successful. However, here is the
> output
> > message:
> > Spectrum decomposition of a symmetric 2nd order tensor
> > *LAPACK error in syevx
> > *Eigenvalues: -0.000236068 0.00423607
> > Spectrum decomposition succeeded!
> >
> > This error message is from "lapack_full_matrix.cc"
> > // Negative return value implies a wrong argument. This should be
> internal.
> > 2049 Assert
> > <
> https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#ga70a0bb353656e704acf927945277bbc6>(info
>
> >= 0, ExcInternalError <
> https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#gab7a0d88175320d08084b2f40f5e3380b
> >());
> > 2050 if (info != 0)
> > 2051 std::cerr << "*LAPACK error in syevx*" << std::endl;
> > **
> > Since I need to run the spectrum decomposition at each Gauss point, it
> is
> > annoying to receive tons of this error message, even though the
> simulation
> > still finishes with no issue.
> >
> > Any insights to avoid this issue?
> >
> > Best,
> >
> > Tao
> >
> > --
> > The deal.II project is located at http://www.dealii.org/
> > <http://www.dealii.org/>
> > For mailing list/f
[deal.II] LAPACK error in syevx
Dear all,
I received the message "LAPACK error in syevx" when using the
*compute_eigenvalues_symmetric* of *lapack_full_matrix* (the interface to
LAPACK in dealii) to compute the eigenvalues of a symmetric matrix with
relatively small matrix entries. I have attached a test file (main.cc) in
the end.
To provide more information, I construct the following 2x2 symmetric matrix
const double *scalar = 1.0;*
for (unsigned int i = 0; i < dim; i++)
for (unsigned int j = 0; j <= i; j++)
{
symmetric_tensor[i][j] = (i + j + 1) * *scalar*;
if (j != i)
symmetric_tensor[j][i] = symmetric_tensor[i][j];
}
*In this case, the eigenvalue decomposition runs without error message:*
Spectrum decomposition of a symmetric 2nd order tensor
Eigenvalues: -0.236068 4.23607
Spectrum decomposition succeeded!
However, if I set
const double *scalar = 1.0e-3; // to make the matrix entries relatively
small*
The eigenvalue decomposition is still successful. However, here is the
output message:
Spectrum decomposition of a symmetric 2nd order tensor
*LAPACK error in syevx*Eigenvalues: -0.000236068 0.00423607
Spectrum decomposition succeeded!
This error message is from "lapack_full_matrix.cc"
// Negative return value implies a wrong argument. This should be internal.
2049 Assert
<https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#ga70a0bb353656e704acf927945277bbc6>(info
>= 0, ExcInternalError
<https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#gab7a0d88175320d08084b2f40f5e3380b>
());
2050 if (info != 0)
2051 std::cerr << "*LAPACK error in syevx*" << std::endl;
Since I need to run the spectrum decomposition at each Gauss point, it is
annoying to receive tons of this error message, even though the simulation
still finishes with no issue.
Any insights to avoid this issue?
Best,
Tao
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/8adfe7d5-220c-4c4a-b3e5-60c2c9c6b278n%40googlegroups.com.
##
# CMake script for the linear thermoelastic monolithic approach:
##
# Set the name of the project and target:
SET(TARGET "main")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0)
FIND_PACKAGE(deal.II 9.4.0
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
/* -
*
* Copyright (C) 2006 - 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* -
*
* Author: Tao Jin
* University of Ottawa, Ottawa, Ontario, Canada
* May 2023
*/
// In this example, the spectrum decomposition of a symmetric 2nd order tensor
// and its relevant derivatives are calculated.
#include
#include
#include
#include
#include
#include
#include
#include
#include
template
void spectrum_decomposition_test()
{
using namespace dealii;
SymmetricT
Re: [deal.II] FullMatrix::copy_to only works for dim <= 3
Hi Wolfgang,
Below is the interface for copy_to:
template
template
void FullMatrix< number >::copy_to ( Tensor< 2, dim > & T,
const size_type src_r_i = 0,
const size_type src_r_j = dim - 1,
const size_type src_c_i = 0,
const size_type src_c_j = dim - 1,
const unsigned int dst_r = 0,
const unsigned int dst_c = 0
) const
When I called the function, I just used the default values for variables 2
to 7.
I actually compared the header file "include/deal.II/lac/full_matrix.h". In
the develop branch (9.5.0), the definition of copy_to is implemented in the
header file. However, in 9.4.0, the header file only has the declaration of
copy_to. The actual implementation happens in "/dealii-9.4.0/src/source/lac/
full_matrix.inst.in"
*for (S : REAL_SCALARS) {template void
FullMatrix::copy_to<1>(Tensor<2, 1> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<2>(Tensor<2, 2> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<3>(Tensor<2, 3> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const; }*
This probably explains why the function works properly in 9.5.0 but does
not work in 9.4.0 for dimensions larger than 3.
Anyhow, I think 9.5.0 has fixed this issue already.
Best,
Tao
On Thursday, June 8, 2023 at 7:36:03 PM UTC-4 Wolfgang Bangerth wrote:
>
> Tao:
>
> On 6/7/23 12:45, Tao Jin wrote:
> > *** Caution: EXTERNAL Sender ***
> >
> > Dear Wolfgang,
> >
> > Here is a simple test code (I also attached the source code):
> >
> > int main()
> > {
> > using namespace dealii;
> >
> > const unsigned int *matrix_dimension* = 4;
> >
> > Tensor<2, matrix_dimension> myTensor;
> >
> > FullMatrix myMatrix(matrix_dimension,
> > matrix_dimension);
> >
> > for (unsigned int i = 0; i < matrix_dimension; i++)
> > for (unsigned int j = 0; j < matrix_dimension; j++)
> > {
> > myMatrix(i,j) = i + j;
> > }
> >
> > myMatrix.copy_to(myTensor);
> >
> > std::cout << myTensor.dimension << std::endl;
> >
> > return 0;
> > }
> >
> > When *matrix_dimension =*2 or 3, the code works fine. When
> *matrix_dimension
> > *= 4, it has a linking error as below:
> > /error: undefined reference to 'void
> > dealii::FullMatrix::copy_to<4>(dealii::Tensor<2, 4, double>&,
> unsigned
> > long, unsigned long, unsigned long, unsigned long, unsigned int,
> unsigned int)
> > const'
> > collect2: error: ld returned 1 exit status
> > make[3]: *** [CMakeFiles/main.dir/build.make:116: main] Error 1/
> >
> > The version of deal.ii is 9.4.0.
>
> I tried this out with the current development version (which will become
> 9.5
> within the next few weeks) and it works just fine:
> https://github.com/dealii/dealii/pull/15331
>
> I am not sure what is going on here. In your code snippet, you call the
> copy_to() function with one argument, but the error message refers to the
> function with 7 arguments. Are you sure that this is the code you wanted
> to use?
>
> Best
> W.
>
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/13481bb3-f449-4972-bb5e-9803460487b9n%40googlegroups.com.
Re: [deal.II] FullMatrix::copy_to only works for dim <= 3
Dear Wolfgang,
Here is a simple test code (I also attached the source code):
int main()
{
using namespace dealii;
const unsigned int *matrix_dimension* = 4;
Tensor<2, matrix_dimension> myTensor;
FullMatrix myMatrix(matrix_dimension,
matrix_dimension);
for (unsigned int i = 0; i < matrix_dimension; i++)
for (unsigned int j = 0; j < matrix_dimension; j++)
{
myMatrix(i,j) = i + j;
}
myMatrix.copy_to(myTensor);
std::cout << myTensor.dimension << std::endl;
return 0;
}
When *matrix_dimension =*2 or 3, the code works fine. When
*matrix_dimension *= 4, it has a linking error as below:
*error: undefined reference to 'void
dealii::FullMatrix::copy_to<4>(dealii::Tensor<2, 4, double>&,
unsigned long, unsigned long, unsigned long, unsigned long, unsigned int,
unsigned int) const'collect2: error: ld returned 1 exit statusmake[3]: ***
[CMakeFiles/main.dir/build.make:116: main] Error 1*
The version of deal.ii is 9.4.0.
I looked into the source code. Indeed, in "full_matrix.h", there is a
declaration of "copy_to". However, the actual implementation is not in
"/dealii-9.4.0/src/source/lac/full_matrix.cc", but in
"/dealii-9.4.0/src/source/lac/full_matrix.inst.in"
*for (S : REAL_SCALARS) {template void
FullMatrix::copy_to<1>(Tensor<2, 1> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<2>(Tensor<2, 2> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<3>(Tensor<2, 3> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const; }*
I assume this is related with explicit instantiation?
Best,
Tao
On Wednesday, June 7, 2023 at 10:38:35 AM UTC-4 Wolfgang Bangerth wrote:
> On 6/6/23 15:46, Tao Jin wrote:
> >
> > However, the above copy_to() function only works when the tensor
> dimension is
> > less or equal to 3, that is,
> > Tensor<2, dim>, dim <=3.
> > When dim is larger than 3, I have to use for loops to copy each matrix
> entry
> > into the rank-2 tensor.
>
> What happens for dim>3? You only say "it only works when", but you don't
> describe what goes wrong if dim>3?
>
>
> > I know that it is typical to set dim as 2 or 3. However, does it make
> sense to
> > make FullMatrix::copy_to() work for a general rank-2 tensor?
>
> I see no reason why this should not work. The function is implemented in
> the
> full_matrix.h file, so I believe it should just work if you call the
> function
> on tensors that are larger than 3x3.
>
> Best
> W.
>
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/447fb768-3783-4dac-805f-89ec14773fc7n%40googlegroups.com.
/* -
*
* Copyright (C) 2006 - 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the licens
[deal.II] FullMatrix::copy_to only works for dim <= 3
Dear all,
When I am implementing a spectrum decomposition for an arbitrary symmetric
matrix, I need to copy a full matrix into a rank-2 tensor. I know that
FullMatrix has a copy_to() member function like below
*templatetemplatevoid FullMatrix< number
>::copy_to ( Tensor< 2, dim > & T,const size_type src_r_i = 0,const
size_type src_r_j = dim - 1,const size_type src_c_i = 0,const size_type
src_c_j = dim - 1,const unsigned int dst_r = 0,const unsigned int dst_c = 0
) const*
However, the above copy_to() function only works when the tensor dimension
is less or equal to 3, that is,
Tensor<2, dim>, dim <=3.
When dim is larger than 3, I have to use for loops to copy each matrix
entry into the rank-2 tensor.
I know that it is typical to set dim as 2 or 3. However, does it make sense
to make FullMatrix::copy_to() work for a general rank-2 tensor?
By the way, in the file full_matrix.inst.in,
*for (S : REAL_SCALARS) {template void
FullMatrix::copy_to<1>(Tensor<2, 1> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<2>(Tensor<2, 2> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const;template void
FullMatrix::copy_to<3>(Tensor<2, 3> &,
const size_type,const
size_type,const size_type,
const size_type,
const unsigned int,
const unsigned int) const; }*
This is probably the reason why copy_to only works for dim <=3
Best,
Tao
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4a697071-4b80-46b7-8f83-3d9df71e3662n%40googlegroups.com.
Re: [deal.II] Form a BlockSparseMatrix based on block matrices from a different BlockSparseMatrix
Dear Wolfgang,
Thank you so much for your guidance. Step-35 is the example where
"SynchronousIterators" is used.
""
using IteratorTuple =
std::tuple::active_cell_iterator,
typename DoFHandler::active_cell_iterator>;
using IteratorPair = SynchronousIterators;
...
WorkStream::run(
IteratorPair(IteratorTuple(dof_handler_velocity.begin_active(),
dof_handler_pressure.begin_active())),
IteratorPair(IteratorTuple(dof_handler_velocity.end(),
dof_handler_pressure.end())),
*this,
::assemble_one_cell_of_gradient,
::copy_gradient_local_to_global,
scratch_data,
per_task_data);
"""
Best,
Tao
On Thursday, April 20, 2023 at 9:22:11 PM UTC-4 Wolfgang Bangerth wrote:
> On 4/20/23 17:27, Tao Jin wrote:
> >
> > My question is, is there a way to involve two different dof_handlers
> when
> > using the WorkStream class, something like the following:
>
> What you are looking for is an object that looks like a *pair* of iterators
> and when you call ++it on it, then it increments both elements of the
> pair.
> The class SynchronousIterator does exactly this.
>
> Alternatively, it is possible to convert an iterator into one DoFHandler
> to an
> iterator into another DoFHandler, like in
> typename DoFHandler::active_cell_iterator cell_stokes
> = stokes_dof_handler.begin_active();
>
> for (..., ++cell_stokes)
> {
> typename DoFHandler::active_cell_iterator
> cell_T (cell_stokes->get_triangulation(),
> cell_stokes->level(),
> cell_stokes->index(),
> _dof_handler);
> ...do something with both iterators...
> }
>
> Best
> W.
>
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/0f1c7da6-773c-46e2-8c3a-f1683085f950n%40googlegroups.com.
Re: [deal.II] Form a BlockSparseMatrix based on block matrices from a different BlockSparseMatrix
Dear Wolfgang,
Thank you so much for your inputs. Indeed, Step-31 is a good example for
what I am looking for. However, I have a follow-up question regarding using
the *WorkStream* class based on TBB during the assembly stage when two
different dof_handlers are involved.
To be specific, in Step-31, two different dof_handlers are defined, one for
the Stoke sub-problem and the other for the thermal sub-problem. Since the
thermal sub-problem depends on the velocity field solved from the Stokes
sub-problem, the following is done to ensure that the cell from the two
dof_handlers are in sync:
auto *cell* = stokes_dof_handler.begin_active();
const auto endc = stokes_dof_handler.end();
auto *temperature_cell* = temperature_dof_handler.begin_active();
for (; cell != endc; ++cell, ++temperature_cell)
{
stokes_fe_values.reinit(*cell*);
temperature_fe_values.reinit(*temperature_cell*);
...
}
My question is, is there a way to involve two different dof_handlers when
using the WorkStream class, something like the following:
WorkStream::run(
m_dof_handler_mechanical.active_cell_iterators(),
[this](const typename DoFHandler::active_cell_iterator ,
ScratchData_ASM_mechanical &
scratch,
PerTaskData_ASM_mechanical &
data) {
this->assemble_system_one_cell_mechanical(cell, scratch, data);
},
[this](const PerTaskData_ASM_mechanical ) {
this->m_constraints_mechanical.distribute_local_to_global(data.m_cell_matrix,
data.m_cell_rhs,
data.m_local_dof_indices,
m_tangent_matrix_mechanical,
m_system_rhs_mechanical);
},
scratch_data,
per_task_data);
Thank you so much!
Best,
Tao
On Sunday, April 16, 2023 at 8:13:41 PM UTC-4 Wolfgang Bangerth wrote:
> On 4/15/23 22:31, Tao Jin wrote:
> >
> > During a time interval, I want to use the first split operator
> > *A B * u1 = b1 *to solve*u1 *and*u2.*
> > *B Cu2b2*
> > Then, plug the solved u1 and u2 into the second split operator
> > *F * u3 = b3*
> > where *F* contains the information of the newly solved *u1* and *u2*.
>
> You will want to look at how step-32 solves the coupled problem. There, it
> is
> a time-dependent one where we first want to solve for the first two
> variables,
> and then the third.
>
> If F depends on u1 and u2, you will have to assemble that matrix separate
> from
> the A,B,C matrices. In that case, you might also want to consider the
> approach
> of step-31.
>
> step-31 and -32 solve the same problem. It might be useful to compare and
> contrast.
>
> Best
> W.
>
> --
>
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/5390f489-a650-40ef-b64d-3689bdf2de72n%40googlegroups.com.
[deal.II] Form a BlockSparseMatrix based on block matrices from a different BlockSparseMatrix
Dear all, This is a similar question related to a post two years ago: https://groups.google.com/g/dealii/c/1L6gLFtvTCo/m/frIyTDXbBwAJ I am developing a fractional time stepping algorithm based on an operator split. Basically, I can obtain the following monolithc system after the finite element discretization: * A B 0u1 b1* *tangent matrix = B C D , solution = u2, RHS = b2* * E 0Fu3 b3* where *tangent matrix* is a 3 by 3 BlockSparseMatrix. During a time interval, I want to use the first split operator *A B * u1 = b1 *to solve* u1 *and* u2.* *B Cu2b2* Then, plug the solved u1 and u2 into the second split operator *F * u3 = b3* where *F* contains the information of the newly solved *u1* and *u2*. My question is, is there an easy way now to extract the sparse matrices *A, B, B, C* to form a new BlockSparseMatrix, something like the following: new BlockSparseMatrix = *tangent matrix*.block(0, 0), *tangent matrix*.block(0, 1) *tangent matrix*.block(1, 0), *tangent matrix*.block(1, 1) Thank you Tao -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/844f4391-91ac-419a-a197-576aa015ab45n%40googlegroups.com.
