Yes, the preconditioner is like a multi-stage pccomposite preconditioner. In one stage, I just precondition a subset of matrix corresponding to certain physical unknowns. To get a better decoupled submatrix, I need to apply that operation.
Thanks. Xiangdong On Wed, Feb 14, 2018 at 9:39 AM, Matthew Knepley <[email protected]> wrote: > On Wed, Feb 14, 2018 at 9:29 AM, Xiangdong <[email protected]> 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. >> > > Barry's point is that this operation is usually rolled into the > preconditioner itself, as in his example of PBJACOBI. Are you building this > preconditioner yourself? > > Matt > > >> Thanks. >> >> Xiangdong >> >> On Tue, Feb 13, 2018 at 6:27 PM, Smith, Barry F. <[email protected]> >> 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 <[email protected]> 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 >>> > >>> > >>> > >>> > >>> >>> >> > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener > > https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/> >
