> On Mar 6, 2018, at 4:37 PM, Ju Liu <liu...@gmail.com> wrote:
> 
> 
> 
> 2018-03-06 14:13 GMT-08:00 Smith, Barry F. <bsm...@mcs.anl.gov>:
> 
> 
> > On Mar 6, 2018, at 4:09 PM, John <johnlucassatur...@gmail.com> wrote:
> >
> > Hi PETSc team:
> >
> > I am trying to create a SIMPLE type approximation for the Schur complement. 
> > It is a nonlinear code, so I have to update the Schur matrix once the 
> > tangent jacobian is updated. I am trying to figure out the best way of 
> > doing that.
> >
> > In terms of generating an algebraic form of Schur S, I am trying to do the 
> > following.
> >
> > // S = A10 x A01, I have performed a diagonal scaling on all matrices, so 
> > the diag(A00) = I.
> > MatMatMult(A10, A01, MAT_INITIAL_MATRIX, S);
> > MatScale(S, -1.0); // S = -1 x S
> > MatAXPY(S, 1.0, A11, S, DIFFERENT_NONZERO_PATTERN); // S = A11 + S
> >
> > My questions are:
> > (1) When I update the Schur complement, do I need to call
> > MatDestroy(&S);
> > MatMatMult(A10, A01, MAT_INITIAL_MATRIX, S);
> > or shall I just call
> > MatMatMult(A10, A01, MAT_REUSE_MATRIX, S);
> 
>    The second form. No reason to delete the matrix and rebuild it each time 
> (delete and rebuild will be a bit slower)
> 
> >
> > Basically, I am not sure how to properly use the flag MAT_REUSE_MATRIX.
> >
> > (2) When I update the Schur complement, the initial call of MatAXPY already 
> > incorporated the contribution of A00 in terms of nonzero patterns. So in my 
> > second call of MatAXPY, shall I use the flag SUBSET_NONZERO_PATTERN?
> 
>    If it is a subset then you should use the subset flag, if they have same 
> nonzero pattern you should use SAME_NONZERO_PATTERN (same nonzero is fastest)
> 
> Thanks Barry! I guess I still need some clarifications on the second point.
> 
> In the first assembly of S, I called MatAXPY to incorporate the entries of 
> A11 into S (with differnet nonzero pattern for safety).
> 
> Now, with new blocks of A, I update S by calling
> MatMatMult(A10, A01, MAT_REUSE_MATRIX, S);

   Ahh, you cannot do this if the MatAXPY introduces new nonzeros into the 
structure of S (the second MatMatMult() will likely crash or generate incorrect 
results since MAT_REUSE_MATRIX requires that each time the routine is called S 
has the exact same nonzero pattern.)

   Now for most problems I would expect that A11 has a subset of the nonzeros 
of S, in that case you can use MatAXPY with subset each time and you can use 
MAT_REUSE_MATRIX. But you would need to verify for your problem that A11 has a 
subset of the nonzeros of S. This is the "ideal" case in terms of performance.

   If A11 does not have a subset of the nonzeros of S then you (obviously) 
cannot use subset and you must destroy the S each time and create a new one 
with each MatMatMult().

   Barry


> At this stage, what is the nonzero pattern of S? If S's pattern is determined 
> by the most recent call of MatMatMult, I guess I will still use the different 
> nonzero pattern. If the pattern of S is inherited from the first assembly, 
> the A11's nonzero pattern is a subset of S's, so subset flag should be safe.
> 
> Thanks!
> 
> John
> 
> 
> 
>  
>    Barry
> 
> >
> > Best regards,
> >
> > John
> 
> 
> 
> 
> -- 
> Ju Liu, Ph.D.
> Postdoctoral Research Fellow
> Stanford School of Medicine
> http://ju-liu.github.io/

Reply via email to