Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Xiangdong]

Hi Barry,

I need some help on the parallel version of MatBAIJBlockDiagonalScale. My
understanding is that MatBAIJBlockDiagonalScale_MPIBAIJ would be a wrapper
on the MatBAIJBlockDiagonalScale_SeqBAIJ. However, I am not clear about how
to get the local part of MPIBAIJ. Does the local part of MPIBAIJ consist of
one or two SeqBAIJ?

Can you show me a similar example of writing a method for MPIBAIJ based on
the SeqBAIJ method? The MatInvertBlockDiagonal is not similar, as that
method only involves the diagonal part A and without the off-diagonal part
B.

Thank you.

Xiangdong

On Wed, Feb 14, 2018 at 2:57 PM, Smith, Barry F.  wrote:

>
>   In the PETSc git branch barry/feature-baij-blockdiagonal-scale I have
> done the "heavy lifting" for what you need. See
> https://bitbucket.org/petsc/petsc/branch/barry/feature-
> baij-blockdiagonal-scale
>
>   It scales the Seq BAIJ matrix by its block diagonal. You will need to
> write a routine to also scale the right hand side vector by the block
> diagonal and then you can try the preconditioner for sequential code. Write
> something like VecBlockDiagonalScale(Vec,const PetscScalar *). You get
> the block size from the vector.
>
>
>   Later you or I can add the parallel version (not much more difficult). I
> don't have time to work on it now.
>
>   Let us know if you have any difficulties.
>
>
>   Barry
>
>
> > On Feb 14, 2018, at 9:10 AM, Xiangdong  wrote:
> >
> > The idea goes back to the alternate-block-factorization (ABF) method
> >
> > https://link.springer.com/article/10.1007/BF01932753
> >
> > and is widely used in the reservoir simulation, where the unknowns are
> pressure and saturation. Although the coupled equations are parabolic, the
> pressure equations/variables are more elliptic and the saturation equations
> are more hyperbolic. People always decouple the transformed linear equation
> to obtain a better (more elliptical) pressure matrix and then apply the AMG
> preconditioner on the decoupled matrix.
> >
> > https://link.springer.com/article/10.1007/s00791-016-0273-3
> >
> > Thanks.
> >
> > Xiangdong
> >
> > On Wed, Feb 14, 2018 at 9:49 AM, Smith, Barry F. 
> wrote:
> >
> >   Hmm, I never had this idea presented to me, I have no way to know if
> it is particularly good or bad. So essentially you transform the matrix
> "decoupling the physics alone the diagonal" and then do PCFIELDSPLIT
> instead of using PCFIELDSPLIT directly on the original equations.
> >
> >   Maybe in the long run this should be an option to PCFIEDLSPLIT. In
> general we like the solvers to manage any transformations, not require
> transformations before calling the solvers. I have to think about this.
> >
> >Barry
> >
> >
> > > On Feb 14, 2018, at 8:29 AM, Xiangdong  wrote:
> > >
> > > The reason for the operation invdiag(A)*A is to have a decoupled
> matrix/physics for preconditioning. For example, after the transformation,
> the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for
> bs=3). One can extract a submatrix (e.g. corresponding to only first
> unknown) and apply special preconditioners for the extracted/decoupled
> matrix. The motivation is that after the transformation, one can get a
> better decoupled matrix to preserve the properties of the unknowns.
> > >
> > > Thanks.
> > >
> > > Xiangdong
> > >
> > > On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F. 
> wrote:
> > >
> > >  In general you probably don't want to do this. Most good
> preconditioners (like AMG) rely on the matrix having the "natural" scaling
> that arises from the discretization and doing a scaling like you describe
> destroys that natural scaling. You can use PCPBJACOBI to use point block
> Jacobi preconditioner on the matrix without needing to do the scaling up
> front. The ILU preconditioners for BAIJ matrices work directly with the
> block structure so again pre-scaling the matrix buys you nothing. PETSc
> doesn't have any particularly efficient routines for computing what you
> desire, the only way to get something truly efficient is to write the code
> directly using the BAIJ data structure, doable but probably not worth it.
> > >
> > >   Barry
> > >
> > >
> > > > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> > > >
> > > > Hello everyone,
> > > >
> > > > I have a block sparse matrices A created from the DMDA3d. Before
> passing the matrix to ksp solver, I want to apply a transformation to this
> matrix: namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the
> block diagonal of A. What is the best way to get the transformed matrix?
> > > >
> > > > At this moment, I created a new mat IDA=inv(diag(A)) by looping
> through each row and call MatMatMult to get B=invdiag(A)*A, then destroy
> the temporary matrix B. However, I prefer the in-place transformation if
> possible, namely, without the additional matrix B for memory saving purpose.
> > > >
> > > > Do 

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Xiangdong]

Thanks a lot, Barry! I see that you had implemented the bs=3 special case.
I will play with these codes and add at least bs=2 case and try to get it
working for parallel baij. I will let you know the update.

Thank you.

Xiangdong

On Wed, Feb 14, 2018 at 2:57 PM, Smith, Barry F.  wrote:

>
>   In the PETSc git branch barry/feature-baij-blockdiagonal-scale I have
> done the "heavy lifting" for what you need. See
> https://bitbucket.org/petsc/petsc/branch/barry/feature-
> baij-blockdiagonal-scale
>
>   It scales the Seq BAIJ matrix by its block diagonal. You will need to
> write a routine to also scale the right hand side vector by the block
> diagonal and then you can try the preconditioner for sequential code. Write
> something like VecBlockDiagonalScale(Vec,const PetscScalar *). You get
> the block size from the vector.
>
>
>   Later you or I can add the parallel version (not much more difficult). I
> don't have time to work on it now.
>
>   Let us know if you have any difficulties.
>
>
>   Barry
>
>
> > On Feb 14, 2018, at 9:10 AM, Xiangdong  wrote:
> >
> > The idea goes back to the alternate-block-factorization (ABF) method
> >
> > https://link.springer.com/article/10.1007/BF01932753
> >
> > and is widely used in the reservoir simulation, where the unknowns are
> pressure and saturation. Although the coupled equations are parabolic, the
> pressure equations/variables are more elliptic and the saturation equations
> are more hyperbolic. People always decouple the transformed linear equation
> to obtain a better (more elliptical) pressure matrix and then apply the AMG
> preconditioner on the decoupled matrix.
> >
> > https://link.springer.com/article/10.1007/s00791-016-0273-3
> >
> > Thanks.
> >
> > Xiangdong
> >
> > On Wed, Feb 14, 2018 at 9:49 AM, Smith, Barry F. 
> wrote:
> >
> >   Hmm, I never had this idea presented to me, I have no way to know if
> it is particularly good or bad. So essentially you transform the matrix
> "decoupling the physics alone the diagonal" and then do PCFIELDSPLIT
> instead of using PCFIELDSPLIT directly on the original equations.
> >
> >   Maybe in the long run this should be an option to PCFIEDLSPLIT. In
> general we like the solvers to manage any transformations, not require
> transformations before calling the solvers. I have to think about this.
> >
> >Barry
> >
> >
> > > On Feb 14, 2018, at 8:29 AM, Xiangdong  wrote:
> > >
> > > The reason for the operation invdiag(A)*A is to have a decoupled
> matrix/physics for preconditioning. For example, after the transformation,
> the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for
> bs=3). One can extract a submatrix (e.g. corresponding to only first
> unknown) and apply special preconditioners for the extracted/decoupled
> matrix. The motivation is that after the transformation, one can get a
> better decoupled matrix to preserve the properties of the unknowns.
> > >
> > > Thanks.
> > >
> > > Xiangdong
> > >
> > > On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F. 
> wrote:
> > >
> > >  In general you probably don't want to do this. Most good
> preconditioners (like AMG) rely on the matrix having the "natural" scaling
> that arises from the discretization and doing a scaling like you describe
> destroys that natural scaling. You can use PCPBJACOBI to use point block
> Jacobi preconditioner on the matrix without needing to do the scaling up
> front. The ILU preconditioners for BAIJ matrices work directly with the
> block structure so again pre-scaling the matrix buys you nothing. PETSc
> doesn't have any particularly efficient routines for computing what you
> desire, the only way to get something truly efficient is to write the code
> directly using the BAIJ data structure, doable but probably not worth it.
> > >
> > >   Barry
> > >
> > >
> > > > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> > > >
> > > > Hello everyone,
> > > >
> > > > I have a block sparse matrices A created from the DMDA3d. Before
> passing the matrix to ksp solver, I want to apply a transformation to this
> matrix: namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the
> block diagonal of A. What is the best way to get the transformed matrix?
> > > >
> > > > At this moment, I created a new mat IDA=inv(diag(A)) by looping
> through each row and call MatMatMult to get B=invdiag(A)*A, then destroy
> the temporary matrix B. However, I prefer the in-place transformation if
> possible, namely, without the additional matrix B for memory saving purpose.
> > > >
> > > > Do you have any suggestion on compute invdiag(A)*A for mpibaij
> matrix?
> > > >
> > > > Thanks for your help.
> > > >
> > > > Best,
> > > > Xiangdong
> > > >
> > > >
> > > >
> > > >
> > >
> > >
> >
> >
>
>

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Smith, Barry F.]

  In the PETSc git branch barry/feature-baij-blockdiagonal-scale I have done 
the "heavy lifting" for what you need. See 
https://bitbucket.org/petsc/petsc/branch/barry/feature-baij-blockdiagonal-scale 

  It scales the Seq BAIJ matrix by its block diagonal. You will need to write a 
routine to also scale the right hand side vector by the block diagonal and then 
you can try the preconditioner for sequential code. Write something like 
VecBlockDiagonalScale(Vec,const PetscScalar *). You get the block size from the 
vector.


  Later you or I can add the parallel version (not much more difficult). I 
don't have time to work on it now. 

  Let us know if you have any difficulties.


  Barry


> On Feb 14, 2018, at 9:10 AM, Xiangdong  wrote:
> 
> The idea goes back to the alternate-block-factorization (ABF) method 
> 
> https://link.springer.com/article/10.1007/BF01932753
> 
> and is widely used in the reservoir simulation, where the unknowns are 
> pressure and saturation. Although the coupled equations are parabolic, the 
> pressure equations/variables are more elliptic and the saturation equations 
> are more hyperbolic. People always decouple the transformed linear equation 
> to obtain a better (more elliptical) pressure matrix and then apply the AMG 
> preconditioner on the decoupled matrix. 
> 
> https://link.springer.com/article/10.1007/s00791-016-0273-3
> 
> Thanks.
> 
> Xiangdong
> 
> On Wed, Feb 14, 2018 at 9:49 AM, Smith, Barry F.  wrote:
> 
>   Hmm, I never had this idea presented to me, I have no way to know if it is 
> particularly good or bad. So essentially you transform the matrix "decoupling 
> the physics alone the diagonal" and then do PCFIELDSPLIT instead of using 
> PCFIELDSPLIT directly on the original equations.
> 
>   Maybe in the long run this should be an option to PCFIEDLSPLIT. In general 
> we like the solvers to manage any transformations, not require 
> transformations before calling the solvers. I have to think about this.
> 
>Barry
> 
> 
> > On Feb 14, 2018, at 8:29 AM, Xiangdong  wrote:
> >
> > The reason for the operation invdiag(A)*A is to have a decoupled 
> > matrix/physics for preconditioning. For example, after the transformation, 
> > the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for 
> > bs=3). One can extract a submatrix (e.g. corresponding to only first 
> > unknown) and apply special preconditioners for the extracted/decoupled 
> > matrix. The motivation is that after the transformation, one can get a 
> > better decoupled matrix to preserve the properties of the unknowns.
> >
> > Thanks.
> >
> > Xiangdong
> >
> > On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F.  wrote:
> >
> >  In general you probably don't want to do this. Most good preconditioners 
> > (like AMG) rely on the matrix having the "natural" scaling that arises from 
> > the discretization and doing a scaling like you describe destroys that 
> > natural scaling. You can use PCPBJACOBI to use point block Jacobi 
> > preconditioner on the matrix without needing to do the scaling up front. 
> > The ILU preconditioners for BAIJ matrices work directly with the block 
> > structure so again pre-scaling the matrix buys you nothing. PETSc doesn't 
> > have any particularly efficient routines for computing what you desire, the 
> > only way to get something truly efficient is to write the code directly 
> > using the BAIJ data structure, doable but probably not worth it.
> >
> >   Barry
> >
> >
> > > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> > >
> > > Hello everyone,
> > >
> > > I have a block sparse matrices A created from the DMDA3d. Before passing 
> > > the matrix to ksp solver, I want to apply a transformation to this 
> > > matrix: namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the 
> > > block diagonal of A. What is the best way to get the transformed matrix?
> > >
> > > At this moment, I created a new mat IDA=inv(diag(A)) by looping through 
> > > each row and call MatMatMult to get B=invdiag(A)*A, then destroy the 
> > > temporary matrix B. However, I prefer the in-place transformation if 
> > > possible, namely, without the additional matrix B for memory saving 
> > > purpose.
> > >
> > > Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?
> > >
> > > Thanks for your help.
> > >
> > > Best,
> > > Xiangdong
> > >
> > >
> > >
> > >
> >
> >
> 
> 

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Xiangdong]

The idea goes back to the alternate-block-factorization (ABF) method

https://link.springer.com/article/10.1007/BF01932753

and is widely used in the reservoir simulation, where the unknowns are
pressure and saturation. Although the coupled equations are parabolic, the
pressure equations/variables are more elliptic and the saturation equations
are more hyperbolic. People always decouple the transformed linear equation
to obtain a better (more elliptical) pressure matrix and then apply the AMG
preconditioner on the decoupled matrix.

https://link.springer.com/article/10.1007/s00791-016-0273-3

Thanks.

Xiangdong

On Wed, Feb 14, 2018 at 9:49 AM, Smith, Barry F.  wrote:

>
>   Hmm, I never had this idea presented to me, I have no way to know if it
> is particularly good or bad. So essentially you transform the matrix
> "decoupling the physics alone the diagonal" and then do PCFIELDSPLIT
> instead of using PCFIELDSPLIT directly on the original equations.
>
>   Maybe in the long run this should be an option to PCFIEDLSPLIT. In
> general we like the solvers to manage any transformations, not require
> transformations before calling the solvers. I have to think about this.
>
>Barry
>
>
> > On Feb 14, 2018, at 8:29 AM, Xiangdong  wrote:
> >
> > The reason for the operation invdiag(A)*A is to have a decoupled
> matrix/physics for preconditioning. For example, after the transformation,
> the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for
> bs=3). One can extract a submatrix (e.g. corresponding to only first
> unknown) and apply special preconditioners for the extracted/decoupled
> matrix. The motivation is that after the transformation, one can get a
> better decoupled matrix to preserve the properties of the unknowns.
> >
> > Thanks.
> >
> > Xiangdong
> >
> > On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F. 
> wrote:
> >
> >  In general you probably don't want to do this. Most good
> preconditioners (like AMG) rely on the matrix having the "natural" scaling
> that arises from the discretization and doing a scaling like you describe
> destroys that natural scaling. You can use PCPBJACOBI to use point block
> Jacobi preconditioner on the matrix without needing to do the scaling up
> front. The ILU preconditioners for BAIJ matrices work directly with the
> block structure so again pre-scaling the matrix buys you nothing. PETSc
> doesn't have any particularly efficient routines for computing what you
> desire, the only way to get something truly efficient is to write the code
> directly using the BAIJ data structure, doable but probably not worth it.
> >
> >   Barry
> >
> >
> > > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> > >
> > > Hello everyone,
> > >
> > > I have a block sparse matrices A created from the DMDA3d. Before
> passing the matrix to ksp solver, I want to apply a transformation to this
> matrix: namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the
> block diagonal of A. What is the best way to get the transformed matrix?
> > >
> > > At this moment, I created a new mat IDA=inv(diag(A)) by looping
> through each row and call MatMatMult to get B=invdiag(A)*A, then destroy
> the temporary matrix B. However, I prefer the in-place transformation if
> possible, namely, without the additional matrix B for memory saving purpose.
> > >
> > > Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?
> > >
> > > Thanks for your help.
> > >
> > > Best,
> > > Xiangdong
> > >
> > >
> > >
> > >
> >
> >
>
>

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Smith, Barry F.]

  Hmm, I never had this idea presented to me, I have no way to know if it is 
particularly good or bad. So essentially you transform the matrix "decoupling 
the physics alone the diagonal" and then do PCFIELDSPLIT instead of using 
PCFIELDSPLIT directly on the original equations. 

  Maybe in the long run this should be an option to PCFIEDLSPLIT. In general we 
like the solvers to manage any transformations, not require transformations 
before calling the solvers. I have to think about this.

   Barry


> On Feb 14, 2018, at 8:29 AM, Xiangdong  wrote:
> 
> The reason for the operation invdiag(A)*A is to have a decoupled 
> matrix/physics for preconditioning. For example, after the transformation, 
> the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for bs=3). 
> One can extract a submatrix (e.g. corresponding to only first unknown) and 
> apply special preconditioners for the extracted/decoupled matrix. The 
> motivation is that after the transformation, one can get a better decoupled 
> matrix to preserve the properties of the unknowns.
> 
> Thanks.
> 
> Xiangdong
> 
> On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F.  wrote:
> 
>  In general you probably don't want to do this. Most good preconditioners 
> (like AMG) rely on the matrix having the "natural" scaling that arises from 
> the discretization and doing a scaling like you describe destroys that 
> natural scaling. You can use PCPBJACOBI to use point block Jacobi 
> preconditioner on the matrix without needing to do the scaling up front. The 
> ILU preconditioners for BAIJ matrices work directly with the block structure 
> so again pre-scaling the matrix buys you nothing. PETSc doesn't have any 
> particularly efficient routines for computing what you desire, the only way 
> to get something truly efficient is to write the code directly using the BAIJ 
> data structure, doable but probably not worth it.
> 
>   Barry
> 
> 
> > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> >
> > Hello everyone,
> >
> > I have a block sparse matrices A created from the DMDA3d. Before passing 
> > the matrix to ksp solver, I want to apply a transformation to this matrix: 
> > namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the block 
> > diagonal of A. What is the best way to get the transformed matrix?
> >
> > At this moment, I created a new mat IDA=inv(diag(A)) by looping through 
> > each row and call MatMatMult to get B=invdiag(A)*A, then destroy the 
> > temporary matrix B. However, I prefer the in-place transformation if 
> > possible, namely, without the additional matrix B for memory saving purpose.
> >
> > Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?
> >
> > Thanks for your help.
> >
> > Best,
> > Xiangdong
> >
> >
> >
> >
> 
> 

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Xiangdong]

The reason for the operation invdiag(A)*A is to have a decoupled
matrix/physics for preconditioning. For example, after the transformation,
the diagonal block is identity matrix ( e.g. [1,0,0;0,1,0;0,0,1]  for
bs=3). One can extract a submatrix (e.g. corresponding to only first
unknown) and apply special preconditioners for the extracted/decoupled
matrix. The motivation is that after the transformation, one can get a
better decoupled matrix to preserve the properties of the unknowns.

Thanks.

Xiangdong

On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F.  wrote:

>
>  In general you probably don't want to do this. Most good preconditioners
> (like AMG) rely on the matrix having the "natural" scaling that arises from
> the discretization and doing a scaling like you describe destroys that
> natural scaling. You can use PCPBJACOBI to use point block Jacobi
> preconditioner on the matrix without needing to do the scaling up front.
> The ILU preconditioners for BAIJ matrices work directly with the block
> structure so again pre-scaling the matrix buys you nothing. PETSc doesn't
> have any particularly efficient routines for computing what you desire, the
> only way to get something truly efficient is to write the code directly
> using the BAIJ data structure, doable but probably not worth it.
>
>   Barry
>
>
> > On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> >
> > Hello everyone,
> >
> > I have a block sparse matrices A created from the DMDA3d. Before passing
> the matrix to ksp solver, I want to apply a transformation to this matrix:
> namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the block
> diagonal of A. What is the best way to get the transformed matrix?
> >
> > At this moment, I created a new mat IDA=inv(diag(A)) by looping through
> each row and call MatMatMult to get B=invdiag(A)*A, then destroy the
> temporary matrix B. However, I prefer the in-place transformation if
> possible, namely, without the additional matrix B for memory saving purpose.
> >
> > Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?
> >
> > Thanks for your help.
> >
> > Best,
> > Xiangdong
> >
> >
> >
> >
>
>

Re: [petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Smith, Barry F.]

 In general you probably don't want to do this. Most good preconditioners (like 
AMG) rely on the matrix having the "natural" scaling that arises from the 
discretization and doing a scaling like you describe destroys that natural 
scaling. You can use PCPBJACOBI to use point block Jacobi preconditioner on the 
matrix without needing to do the scaling up front. The ILU preconditioners for 
BAIJ matrices work directly with the block structure so again pre-scaling the 
matrix buys you nothing. PETSc doesn't have any particularly efficient routines 
for computing what you desire, the only way to get something truly efficient is 
to write the code directly using the BAIJ data structure, doable but probably 
not worth it.

  Barry


> On Feb 13, 2018, at 5:21 PM, Xiangdong  wrote:
> 
> Hello everyone,
> 
> I have a block sparse matrices A created from the DMDA3d. Before passing the 
> matrix to ksp solver, I want to apply a transformation to this matrix: namely 
> A:= invdiag(A)*A. Here invdiag(A) is the inverse of the block diagonal of A. 
> What is the best way to get the transformed matrix?
> 
> At this moment, I created a new mat IDA=inv(diag(A)) by looping through each 
> row and call MatMatMult to get B=invdiag(A)*A, then destroy the temporary 
> matrix B. However, I prefer the in-place transformation if possible, namely, 
> without the additional matrix B for memory saving purpose.
> 
> Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?
> 
> Thanks for your help.
> 
> Best,
> Xiangdong
> 
> 
> 
>  

[petsc-users] multiply a mpibaij matrix by its block diagonal inverse

[Xiangdong]

Hello everyone,

I have a block sparse matrices A created from the DMDA3d. Before passing
the matrix to ksp solver, I want to apply a transformation to this matrix:
namely A:= invdiag(A)*A. Here invdiag(A) is the inverse of the block
diagonal of A. What is the best way to get the transformed matrix?

At this moment, I created a new mat IDA=inv(diag(A)) by looping through
each row and call MatMatMult to get B=invdiag(A)*A, then destroy the
temporary matrix B. However, I prefer the in-place transformation if
possible, namely, without the additional matrix B for memory saving purpose.

Do you have any suggestion on compute invdiag(A)*A for mpibaij matrix?

Thanks for your help.

Best,
Xiangdong