[petsc-users] Copy MPIAIJ Matrix onto MPISeq Matrices on All Processors

Tue Mar 20 10:06:34 CDT 2018

    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 <aliberkkahraman at yahoo.com> 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 <aliberkkahraman at yahoo.com> 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 <aliberkkahraman at yahoo.com> 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 <aliberkkahraman at yahoo.com> 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/
>