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 <epsco...@gmail.com> 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. <bsm...@mcs.anl.gov> 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 <epsco...@gmail.com> 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
> >
> >
> >
> >
> 
> 

Reply via email to