I made the change in the function so I now get this error PetscErrorCode ComputeRho(DMMG *dmmg, Vec rho) { DA distA = (DA)dmmg->dm; In function ‘PetscErrorCode Myfunc(_n_DMMG**, _p_Vec*)’: error: request for member ‘dm’ in ‘* dmmg’, which is of non-class type ‘_n_DMMG*’ It looks like I have a null pointer to _n_DMMG** instead of _n_DMMG* Chris On Mon, Nov 22, 2010 at 1:34 PM, Barry Smith <[email protected]> wrote:

On Nov 22, 2010, at 3:23 PM, Chris Bording wrote:

...

Hi I am trying to pass the DMMG struct to a function similar to what is in the examples.

extern PetscErrorCode Myfunc(DMMG,Vec);

int main(int argc,char **argv) { DMMG *dmmg; ^^^^^^ Note this ...Vec b;

Myfunc(dmmg,b) }

PetscErrorCode Myfunc(DMMG *dmmg, Vec b) ^^^ add this

...

{ DA da = (DA)dmmg->dm; ..... }

error: cannot convert ‘_n_DMMG**’ to ‘_n_DMMG*’ for argument ‘1’ to ‘PetscErrorCode Myfunc(_n_DMMG*, _p_Vec*)

So how do I pass a DMMG defined structure into my function?

Thanks Chris

On Tue, Nov 23, 2010 at 00:27, Sean Farley <[email protected]> wrote:

2) Change your function declaration back to what you had:

PetscErrorCode ComputeRho(DMMG dmmg, Vec rho) {...}

and then pass the object to your function as so:

ComputeRho(*dmmg, b);

Do you need your dmmg to be a pointer in your 'main' function? It looks like your declaration could just be:

DMMG is simultaneously an array representing a full hierarchy and the individual levels in the hierarchy. If you have DMMG *dmmg; then dmmg[0] (same as *dmmg) represents the coarsest level. If your function only works on one level, then you should just pass in a DMMG (not a DMMG*). If your function only needs the dmmg->dm field, then just make your function use a DA, then make that the argument of your function and call it as DMMG *dmmg; // dmmg is initialized somehow ... MyFunc(DMMGGetDA(dmmg),bvec); This passes in the DA on the finest level. DMMG is a weird beast. Once its functionality has been migrated to SNES and KSP (underway in petsc-dev), it will be removed. Jed Jed

On Tue, Nov 23, 2010 at 10:32, Clemens Domanig <[email protected]>wrote:

I'm writing a FEM-program and at the moment I try to parallize the assembly of the stiffness-matrix. The idea would be if I for example use 4 threads I create 4 matrices and 4 rhs-vectors. Each thread fills its matrix and rhs-vector. At the end I add all matrices and vectors.

Adding sparse matrices together is a relatively expensive operation, usually more expensive than assembling the matrix you wanted in the first place.

But the Petsc-manual says that Petsc is not thread-safe.

There is a minor issue of some logging functions that are not thread safe. Wrapping locks around some of those operations would not be a huge deal, but hasn't been done.

Is that the reason why it doesn't work?

The matrix data structures cannot cheaply be mutated in a thread-safe way. We would either need fine-grained locks which have more overhead, and are very tricky when the user does not preallocate correctly (requiring "rollback" logic), or a coarse-grained lock for the whole MatSetValues. If you want to do multi-threaded assembly, you should just put your own lock around your call to MatSetValues. As long as your physics does some work (e.g. is a bit more than a linear constant-coefficient problem on affine elements), this can work fine for a few threads, usually 8-16 or so, before serialization of MatSetValues becomes a bottleneck. Note that you can also use multiple MPI processes to keep your cores busy, there have been a few discussions on this list, and petsc-dev, about the relative merits of threads versus processes. If you tell us a bit about your problem, we may be able to predict how each will perform. This helps to assess the importance of making PETSc thread-safe and making certain kernels perform well with threads, relative to other features. Jed