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 <[email protected]> wrote: > On Sun, Apr 27, 2014 at 9:05 AM, Justin Dong <[email protected]> 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 <[email protected]>wrote: >> >>> On Sun, Apr 27, 2014 at 7:47 AM, Justin Dong <[email protected]> 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 <[email protected]>wrote: >>>> >>>>> On Sun, Apr 27, 2014 at 6:25 AM, Justin Dong <[email protected]> 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 >>>>>> <[email protected]>wrote: >>>>>> >>>>>>> On Sun, Apr 27, 2014 at 4:14 AM, Justin Dong <[email protected]> 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 >
