2018-03-06 14:49 GMT-08:00 Smith, Barry F. <bsm...@mcs.anl.gov>: > > > > 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(). > > I see. Thanks a lot!
John > 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/ > > -- Ju Liu, Ph.D. Postdoctoral Research Fellow Stanford School of Medicine http://ju-liu.github.io/