[petsc-users] Assembling a finite element matrix in parallel
Justin Dong
jsd1 at rice.edu
Sun Apr 27 19:22:33 CDT 2014
Thanks for your help -- the matrix assembly is working correctly now.
However, say I want to solve the linear problem in serial now. I have
A_global and b_global created using MatCreateAIJ and VecCreateMPI, but want
to solve A_global*x = b_global using PCLU. I have:
KSPCreate(PETSC_COMM_SELF, &ksp);
KSPSetOperators(ksp, A_global, A_global, DIFFERENT_NONZERO_PATTERN);
KSPGetPC(ksp, &pc);
PCSetType(pc, PCLU);
KSPSetTolerances(ksp, 1.e-13, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
KSPSolve(ksp, b_global, b_global);
Is there anyway to get this to work? I'd prefer using PCLU rightnow, and
besides, the speed-up from assembling in parallel is already significant.
On Sun, Apr 27, 2014 at 9:11 AM, Matthew Knepley <knepley at gmail.com> wrote:
> On Sun, Apr 27, 2014 at 9:05 AM, Justin Dong <jsd1 at rice.edu> wrote:
>
>> When I try to invoke MatCreateMPIAIJ, I'm getting this warning:
>>
>> *main.c:1842:2: **warning: **implicit declaration of function
>> 'MatCreateMPIAIJ' is invalid in C99 [-Wimplicit-function-declaration]*
>>
>> MatCreateMPIAIJ(PETSC_COMM_SELF, NLoc, NLoc, NLoc*nElems,
>> NLoc*nElems, (ref+3)*NLoc,
>>
>> * ^*
>>
>> and this error:
>>
> The name was changed in the last release:
>
>
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateAIJ.html
>
> Thanks,
>
> Matt
>
>
>> Undefined symbols for architecture x86_64:
>>
>> "_MatCreateMPIAIJ", referenced from:
>>
>> _DG_SolveOil in main.o
>>
>> ld: symbol(s) not found for architecture x86_64
>>
>> clang: error: linker command failed with exit code 1 (use -v to see
>> invocation)
>>
>> make: [main] Error 1 (ignored)
>>
>> /bin/rm -f -f main.o
>>
>>
>> This is how I'm calling MatCreateMPIAIJ:
>>
>> MatCreateMPIAIJ(PETSC_COMM_WORLD, nlocal, nlocal, nlocal*nelements,
>> nlocal*nelements, 3*nlocal,
>>
>> PETSC_NULL, 3*nlocal, PETSC_NULL, &A_global);
>>
>>
>> I'm able to run some simple MPI functions and use mpiexec, so I don't
>> think it's an issue with MPI itself.
>> Sincerely,
>> Justin
>>
>>
>>
>> On Sun, Apr 27, 2014 at 8:35 AM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Sun, Apr 27, 2014 at 7:47 AM, Justin Dong <jsd1 at rice.edu> wrote:
>>>
>>>> I think I've got the right idea on how to divide up the for loop:
>>>>
>>>> MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
>>>> MPI_Comm_rank(MPI_COMM_WORLD,&myid);
>>>>
>>>> mystart = (nelements / numprocs) * myid;
>>>> if (nelements % numprocs > myid)
>>>> {
>>>> mystart += myid;
>>>> myend = mystart + (nelements / numprocs) + 1;
>>>> }
>>>> else
>>>> {
>>>> mystart += nelements % numprocs;
>>>> myend = mystart + (nelements / numprocs);
>>>> }
>>>>
>>>
>>> We have a function that does exactly this:
>>> http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Sys/PetscSplitOwnership.html
>>>
>>>
>>>> But then, do I do this?
>>>>
>>>
>>> Yes.
>>>
>>>
>>>> for (k = mystart+1; k < myend; ++k)
>>>> {
>>>> get_index(k,ie,je); /* ie and je are arrays that indicate where to
>>>> place A_local */
>>>> compute_local_matrix(A_local,...);
>>>>
>>>> MatSetValues(A_global, nlocal, ie, nlocal, je, A_local, ADD_VALUES);
>>>> }
>>>>
>>>> The indices are global indices and I'm doing
>>>>
>>>> MatCreateSeqAIJ(PETSC_COMM_WORLD, nglobal, nglobal, 3*nlocal,
>>>> PETSC_NULL, &A_global);
>>>>
>>>
>>> This should be MatCreateMPIAIJ() for parallel execution.
>>>
>>> Matt
>>>
>>> to create the matrix. Running the program seems to give me multiple
>>>> errors, but mainly
>>>>
>>>> [3]PETSC ERROR: Object is in wrong state!
>>>>
>>>> [3]PETSC ERROR: Mat object's type is not set: Argument # 1!
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Sun, Apr 27, 2014 at 6:54 AM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>
>>>>> On Sun, Apr 27, 2014 at 6:25 AM, Justin Dong <jsd1 at rice.edu> wrote:
>>>>>
>>>>>> Hi Matt,
>>>>>>
>>>>>> For option 1), that would be using MPI and not any special functions
>>>>>> in PETSc? I ask since I've never done parallel programming before. I tried
>>>>>> using OpenMP to parallelize that for loop but it resulted in some errors
>>>>>> and I didn't investigate further, but I'm assuming it's because each
>>>>>> process wasn't communicating properly with the others during MatSetValues?
>>>>>>
>>>>>
>>>>> Yes, using MPI, so each process owns a set of elements that it loops
>>>>> over. The Mat object manages the global
>>>>> values as long as you use global indices for the (row, column). There
>>>>> are of course many refinements for this,
>>>>> but I think the thing to do is get something working fast, and then
>>>>> optimize it.
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>> Sincerely,
>>>>>> Justin
>>>>>>
>>>>>>
>>>>>> On Sun, Apr 27, 2014 at 5:57 AM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>>
>>>>>>> On Sun, Apr 27, 2014 at 4:14 AM, Justin Dong <jsd1 at rice.edu> wrote:
>>>>>>>
>>>>>>>> Hi all,
>>>>>>>>
>>>>>>>> I'm currently coding a finite element problem in PETSc. I'm
>>>>>>>> computing all of the matrices by myself and would prefer to do it that way.
>>>>>>>> I want to parallelize the assembly of the global matrices. This is a
>>>>>>>> simplified version of the assembly routine (pseudocode):
>>>>>>>>
>>>>>>>> for (k = 0; k < nelements; ++k)
>>>>>>>> {
>>>>>>>> get_index(k,ie,je); /* ie and je are arrays that indicate where
>>>>>>>> to place A_local */
>>>>>>>> compute_local_matrix(A_local,...);
>>>>>>>>
>>>>>>>> MatSetValues(A_global, nlocal, ie, nlocal, je, A_local,
>>>>>>>> ADD_VALUES);
>>>>>>>> }
>>>>>>>>
>>>>>>>> This is for DG finite elements and I know the matrix easily be
>>>>>>>> assembled in parallel. Even on my laptop, this would allow me to
>>>>>>>> significantly speed up my code. The only problem is, I'm not sure how to do
>>>>>>>> this in PETSc. I'd assume I need to use MatCreateMPIAIJ instead of
>>>>>>>> MatCreateSeqAIJ, but wasn't able to find any using tutorials on how I might
>>>>>>>> do this.
>>>>>>>>
>>>>>>>
>>>>>>> 1) If you just split this loop across processes, it would work
>>>>>>> immediately. However, that can be non-optimal in terms
>>>>>>> of communication.
>>>>>>>
>>>>>>> 2) PETSc matrices are distributed by contiguous blocks of rows (see
>>>>>>> manual section on matrices), so you would like
>>>>>>> those row blocks to correspond roughly to your element blocks.
>>>>>>>
>>>>>>> 3) You will also have to preallocate, but that should be the same as
>>>>>>> you do now for the sequential case, except you
>>>>>>> check whether a column is inside the diagonal block.
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> Matt
>>>>>>>
>>>>>>>
>>>>>>>> Sincerely,
>>>>>>>> Justin
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> 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
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> 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
>>>
>>
>>
>
>
> --
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140427/e29965b1/attachment.html>
More information about the petsc-users
mailing list