Thank you for all the help. I think I'm mostly clear now. I think my only remaining question relates to point 1: So let's call the global matrix across all processors A and let's say it's a square matrix of size n. Then each processor has a subset of the rows of A so that there's in principle a contiguous rectangular chunk of A on each proc. Now we want to use this preallocator matrix to determine the non-zero structure of each processors chunk of A. Let's call this matrix P. You say that P should also be square. But I would think that, at least locally, P should be rectangular on each proc and basically have the same size as the local chunks of A. Do you mean that P is square when considered globally across all processors? I guess my question is a bit subtle, but basically, is P essentially a parallel Mat type which exists nominally across all procs, or is it a serial type object which is local to each proc? I think I was viewing it in the second way since it doesn't obviously require communication but I guess if you consider it in the first way it makes sense that it would be nxn.
On Fri, Apr 1, 2022 at 12:00 PM Matthew Knepley <knep...@gmail.com> wrote: > On Fri, Apr 1, 2022 at 12:57 PM Samuel Estes <samueleste...@gmail.com> > wrote: > >> Ah, so let me see if I understand this now. So basically you preallocate >> for this "dummy matrix" of type MatPreallocator. But rather than simply >> calling MatXXXAIJSetPreallocation() to allocate space, you actually tell >> this "dummy matrix" exactly where the non-zeros will be. You then call this >> MatPreallocatorPreallocate() routine to pass the non-zero structure of the >> "dummy matrix" to the actual matrix that you're using for the linear >> system. You then call MatSetValues() (or a variant thereof) to set the >> values of this preallocated matrix. >> >> A couple more clarifying questions: >> 1. So this MatPreallocator matrix should have a rectangular structure of >> (rEnd-rStart) rows and n columns (where n is the total size of the system? >> I'm assuming here that the global matrix across all processors is square. >> > > This "dummy matrix" will look exactly like your system matrix, meaning it > is square, and you tell it the nonzeros by just calling MatSetValues(), > exactly as you would for your system matrix. That is why you can reuse > exactly the same code. > > >> 2. You preallocate for the dummy matrix by calling matsetvalues rather >> than the usual way of preallocating? >> > > Yes. > > >> 3. That's what you mean by looping twice? You loop over once using >> matsetvalues to preallocate, then you have to loop using matsetvalues again >> a second time to actually set the values of the parallel Mat you actually >> use to solve the system? >> > > Yes. > > Thanks, > > Matt > > >> On Fri, Apr 1, 2022 at 11:50 AM Matthew Knepley <knep...@gmail.com> >> wrote: >> >>> On Fri, Apr 1, 2022 at 12:45 PM Samuel Estes <samueleste...@gmail.com> >>> wrote: >>> >>>> Thanks! This seems like it might be what I need. I'm still a little >>>> unclear on how it works though? My problem is basically that for any given >>>> row, I know the total non-zeros but not how many occur in the diagonal vs >>>> off-diagonal. Without knowledge of the underlying grid, I'm not sure how >>>> there could be a black box utility to figure this out? Am I >>>> misunderstanding how this is used? >>>> >>> >>> So each process gets a stack of rows, [rStart, rEnd). A nonzero (r, c) >>> is in the diagonal block if rStart <= c < rEnd. So if you know (r, c) for >>> each nonzero, you know whether it is in the diagonal block. >>> >>> Thanks, >>> >>> Matt >>> >>> >>>> On Fri, Apr 1, 2022 at 11:34 AM Matthew Knepley <knep...@gmail.com> >>>> wrote: >>>> >>>>> On Fri, Apr 1, 2022 at 12:27 PM Samuel Estes <samueleste...@gmail.com> >>>>> wrote: >>>>> >>>>>> Hi, >>>>>> >>>>>> I have a problem in which I know (roughly) the number of non-zero >>>>>> entries for each row of a matrix but I don't have a convenient way of >>>>>> determining whether they belong to the diagonal or off-diagonal part of >>>>>> the >>>>>> parallel matrix. Is there some way I can just allocate the total number >>>>>> of >>>>>> non-zeros in a row regardless of which part they belong to? I'm assuming >>>>>> that this is not possible but I just wanted to check. It seems like it >>>>>> should be possible in principle since the matrix is only split by the >>>>>> rows >>>>>> but the columns of a row are all on the same processor (at least as I >>>>>> understand it). Thanks! >>>>>> >>>>> >>>>> In serial, the matrix is stored by rows. In parallel, it is split into >>>>> a diagonal and off-diagonal block, so that we can overlap communication >>>>> and >>>>> computation in the matvec. >>>>> >>>>> However, we have a convenience structure for figuring this out, called >>>>> MatPreallocator, >>>>> https://petsc.org/main/docs/manualpages/Mat/MATPREALLOCATOR.html >>>>> In my code, I wrote a loop around the code that filled up my matrix, >>>>> which executed twice. On the first pass, I fed in the MatPreallocator >>>>> matrix. When this finished >>>>> you can call >>>>> https://petsc.org/main/docs/manualpages/Mat/MatPreallocatorPreallocate.html#MatPreallocatorPreallocate >>>>> on your system amtrix, then on the second >>>>> pass use the system matrix. This was only a few extra lines of code >>>>> for me. If you want to optimize further, you can have a flag that only >>>>> computes the values on the >>>>> second pass. >>>>> >>>>> Thanks, >>>>> >>>>> Matt >>>>> >>>>> Sam >>>>>> >>>>> -- >>>>> 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.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/ >>> <http://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/ > <http://www.cse.buffalo.edu/~knepley/> >