The FD construction of Jacobians with coloring in PETSc requires that you
provide the nonzero structure of the Jacobian; that is you need to provide the
Jacobian but without the numerical values. PETSc coloring routines can then
color the matrix and then PETSc can compute the Jacobian numerical values.
Here is a chunk of code that does what you need
/* Create coloring context */
ierr = MatColoringCreate(Jac,&mc);CHKERRQ(ierr);
ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr);
ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
ierr = MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);CHKERRQ(ierr);
ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode
(*)(void))FormFunction,&user);CHKERRQ(ierr);
ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
ierr = MatFDColoringSetUp(Jac,iscoloring,matfdcoloring);CHKERRQ(ierr);
/* ierr =
MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
ierr =
SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr);
ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
Barry
> On Mar 20, 2018, at 9:56 AM, Ali Berk Kahraman <[email protected]>
> wrote:
>
> Now I see my problem. What I want to do is to generate many rows of Jacobian
> once. For this, the method that must be used is coloring. If I am to try it
> without calling coloring, I have no way but to give the all processes the
> whole data, but as you have pointed out it is very very nonscalable. It
> instantly multiplies the memory use with number of processors.
>
> Thank you very much for your constructive replies.
>
> Ali
>
> On 20-03-2018 17:30, Matthew Knepley wrote:
>> On Tue, Mar 20, 2018 at 10:11 AM, Ali Berk Kahraman
>> <[email protected]> wrote:
>> First of all, thank you for elaboration.
>>
>> The problem is, I want to generate J myself, using FD, for the following
>> reason.
>> Okay. This generation should happen entirely inside FormJacobian(), so you
>> do not need MFFD.
>>
>> I have a WORLD matrix M that I feed into F(u) using user defined context
>> pointer *ctx. And I call collective functions on that matrix.
>> That should be fine. To be very precise, when you create the SNES, you give
>> it a comm. The FormFunction() and FormJacobian()
>> functions are collective on this comm.
>>
>> I call a for loop inside CreateJacobian function to perturb the F() entries
>> one by one. But since I have collective calls on a WORLD matrix M inside F,
>> I have to do this virtually sequentially(not "truly" sequential, but other
>> processes must wait while one is calling F()). I cannot have two different
>> processes perturb the u vector at different locations, this confuses the
>> collective calls on that matrix (or so I have observed, maybe there was
>> another bug in my code.).
>>
>> Okay, this is a "bug" in your FD algorithm. What collective calls are you
>> making when doing the finite difference?
>>
>> Note that if you really do need collective calls, then you are serializing
>> the construction of this matrix, meaning no matter
>> how many procs you use, the sequence will still be as long as one proc. If,
>> however, that is really what you want, then
>> you need to make sure that everyone calls that function for each
>> perturbation. This means you need to tell everyone how
>> many calls they need to make (either sum(perturbations) or
>> max(perturbations) depending on what you are doing). You
>> can use MPI_Allreduce() for this at the start of your function.
>> Here, what I want to do is to have a SELF copy of the u vector and the input
>> matrix M to F() on every process, so every process can independently work
>> the data. So, I need a way to copy the WORLD matrix into SELF. My current
>> idea is to write the WORLD matrix as a binary file, and read it into a SELF
>> matrix.
>> That is really nonscalable. I cannot believe this is really what you want to
>> do. However, you can do it using
>> VecScatterCreateToAll() and
>>
>>
>> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateSubMatrices.html
>>
>> However, I would write out the algorithm in more detail than you have
>> explained it here. So far, I am not
>> convinced you should be doing this.
>>
>> Thanks,
>>
>> Matt
>>
>> I am aware that what I want to do should be done using coloring, much memory
>> and time effective. But, I am working on adaptive grid finite differences so
>> there is not an easy way to create coloring on my grid. If there is a Petsc
>> function that generates the coloring for me, I can also try it.
>> My apologies, even the simple form got this long :)
>>
>>
>> On 20-03-2018 16:47, Matthew Knepley wrote:
>>> On Tue, Mar 20, 2018 at 9:31 AM, Ali Berk Kahraman
>>> <[email protected]> wrote:
>>> Thank your for the reply.
>>>
>>> I have tried it, but could not get it to work. For it to work, I have to
>>> know at which vector entry the epsilon in the FD approximation is to be
>>> added. It is an input for my inexpensive function, and I could not find a
>>> way to pass it in using MFFD context.
>>>
>>> MFFD has it in the context as current_u. There is no accessor, but we can
>>> easily write one. However, that is not how MFFD works. Here is the idea:
>>>
>>> You give us a function F, your inexpensive function. We get to call F
>>> however we want to call it, and produce the action of the matrix. So for
>>> instance,
>>> we could do this
>>>
>>> F(x + h u) - F(x)
>>> ---------------------
>>> h ||u||
>>>
>>> but we do not have to do that. You should not worry about how we perturb
>>> the function, you just return the value of the function when we call it.
>>>
>>> I have been trying to get it working with PETSC_COMM_WORLD too with no
>>> luck. My guess is that my FormRHSFunction becomes a collective function on
>>> WORLD, so I cannot generate many separate copies each doing its own thing.
>>> For this reason I am trying to move it to SELF, but I cannot do that if I
>>> cannot find a way to represent my WORLD matrix input fully on each process.
>>>
>>> I think things have gotten too confused here for us to help. Lets step back
>>> and simplify. Can you tell us briefly exactly what you want to accomplish?
>>> I will start by
>>> explaining SNES. We are essentially doing Newton's method. This means
>>> calculating a step du, based on an initial solution u, using this equation
>>>
>>> J(u) du = -F(u)
>>>
>>> so we need a function to form the residual F, and maybe the Jacobian J but
>>> we can also form it using calls to F.
>>>
>>> It sounds like you want to do something inside F, but I cannot figure out
>>> what yet.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>> Sorry I am stuck in the same question, but any ideas on how to move a
>>> matrix from WORLD to SELF, similar to the way it is possible with vectors
>>> using VecScatter?
>>>
>>>
>>>
>>> On 19-03-2018 21:08, Smith, Barry F. wrote:
>>>
>>> On Mar 19, 2018, at 12:00 PM, Ali Berk Kahraman <[email protected]>
>>> wrote:
>>>
>>> Dear All,
>>>
>>> I have to create my Jacobian numerically using an FD function I have
>>> written, but the problem is that the function evaluation FormRHSFunction
>>> creates PETSC_COMM_WORLD objects in it, and takes 1 PETSC_COMM_WORLD matrix
>>> as an input. So I cannot work it on parallel (every processor calculating
>>> its own portion of Jacobian), I tried but the machine gets confused.
>>> Huh. It is fine to have the function work on PETSC_COMM_WORLD, and it
>>> can be in parallel.
>>> My initial idea was to create the objects in PETSC_COMM_SELF type within
>>> the RHS function, but the input COMM_WORLD matrix standed where it was. So
>>> my question is, is there a way to get a complete copy of a distributed
>>> matrix on all the processors as seq matrices?
>>>
>>> I cannot use the Petsc's version FD function because the FormRHSFunction I
>>> call to create the Jacobian is slightly different than the original
>>> FormRHSFunction. The original function is too expensive, so I only
>>> calculate the relevant parts of it in the Jacobian function.
>>> You can do this with MatCreateMFFD() and set your inexpensive function
>>> as the function.
>>>
>>> Best Regards,
>>>
>>> Ali Berk Kahraman
>>> M.Sc. Student, Mechanical Engineering
>>> Bogazici Uni. Istanbul, Turkey
>>>
>>>
>>>
>>>
>>>
>>> --
>>> 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/
>>
>>
>>
>>
>> --
>> 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/
>
