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 > > > > > > > > > >
