libMesh: A C++ Finite Element Library

Dear Ben,

On Fri, 17 Apr 2009, Kirk, Benjamin (JSC-EG311) wrote:

> Please find attached a patch which addresses the ridiculous implementation
> of adding a scalar to a PetscVector pointed out last week...  Let me know if
> it works for you, if so I'll submit it.

Actually, I'm not using that function anywhere in my code, so there is 
no easy way for me to test it.  (Remember that Jed was the one who 
pointed this out initially.)

On the other hand, there are some more possibly inefficient things in 
that file.  For instance, PetscVector::insert() calls 
PetscVector::set() for each index, and I suppose it would be better to 
use VecSetValues() instead.  It's difficult to decide how far one 
should go at the moment with optimizing the PetscVector class.

Anyway, I looked over your patch and it seems correct to me, and in 
particular it coincides with Jed's suggestion.  Since Jed seems to be 
familiar with PETSc very well, this could be considered of being 
enough reason to submit the patch right now.

Other optimizations should probably wait until we know whether this is 
really a bottleneck.  Before I run the application with the PETSc log 
output option that Jed suggested, I would like the remining bug to be 
fixed, since one never knows what that implies.  I'm currently waiting 
for Roy to report whether he can reproduce the bug; I guess that he 
has been too busy by now.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Fri, 17 Apr 2009, Tim Kroeger wrote:

> I would like the remining bug to be fixed, since one never knows
> what that implies.  I'm currently waiting for Roy to report whether
> he can reproduce the bug; I guess that he has been too busy by now.

Sorry about the delay; you guess right.

The prodding was helpful.  I did have time to set up your test program
and get it started this morning.  Not sure if/when it's going to
finish, though.  I'd like to run all 24 processes on the same box for
ease of debugging, but it's currently halfway through the refinement
steps and using 6GB memory out of 4GB available RAM.  We've got a
couple 16GB nodes here, but they're even busier than I am through most
of the week lately.  But if the 4GB server proves insufficient I'll
see if I can monopolize one of the big guys tomorrow or Sunday.
---
Roy

Dear Roy,

On Fri, 17 Apr 2009, Roy Stogner wrote:

> I did have time to set up your test program and get it started this 
> morning.  Not sure if/when it's going to finish, though.

Did it reveal anything yet?

> I'd like to run all 24 processes on the same box for
> ease of debugging, but it's currently halfway through the refinement
> steps and using 6GB memory out of 4GB available RAM.

You mean, it was already swapping?  I see, that's going to slow it 
down essentially.

By the way: Nobody seems to have checked in my patch that I sent to 
the list last week (April 15) (nor has anybody stated that my patch 
could cause problems).

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Wed, 22 Apr 2009, Tim Kroeger wrote:

> On Fri, 17 Apr 2009, Roy Stogner wrote:
>
>> I did have time to set up your test program and get it started this 
>> morning.  Not sure if/when it's going to finish, though.
>
> Did it reveal anything yet?

It left the system hosed for long enough that I had to kill it.  I'll
be able to try again on a bigger system this weekend.

> By the way: Nobody seems to have checked in my patch that I sent to the list 
> last week (April 15) (nor has anybody stated that my patch could cause 
> problems).

It's one of the many things backing up in my inbox right now, sorry.
---
Roy

On Wed, 22 Apr 2009, Roy Stogner wrote:

> On Wed, 22 Apr 2009, Tim Kroeger wrote:
>
>> By the way: Nobody seems to have checked in my patch that I sent to the 
>> list last week (April 15) (nor has anybody stated that my patch could cause 
>> problems).
>
> It's one of the many things backing up in my inbox right now, sorry.

Patch looks good.  I'm in a rush now, but I'll add it late tonight or
tomorrow.  One question first: can anyone think of a better name than
constrain_nothing()?  The analogy to our other constrain_* functions
makes sense, but it seems odd outside that context to have a method
essentially named "do_nothing" that does modify its inputs.

I can't think of any better name myself, so I'll commit the new
method as constrain_nothing().  But if someone *can* think of a more
descriptive name, let me know so we can change it while the API's
still only got one user.  ;-)
---
Roy

Dear Roy,

On Wed, 22 Apr 2009, Roy Stogner wrote:

> One question first: can anyone think of a better name than 
> constrain_nothing()?  The analogy to our other constrain_* functions 
> makes sense, but it seems odd outside that context to have a method 
> essentially named "do_nothing" that does modify its inputs.

Well, of course you are right.  But actually, I find this is only a 
symptom of the fact that a user might not expect the constrain_*() 
methods to modify their dof_indices argument.  At least, that was true 
for me quite a long time, and it caused a number of programming errors 
in my applications.  I learned this when I implemented the 
constrain_dyad_matrix() function.  Since that time I understand that 
these methods *have* to do this.

What I want to say is this: First, having a method called 
"constrain_nothing()" might make the innocent user look into the 
details earlier and prevent him from making mistakes.  Second, the 
whole constaining mechanism could be reworked completely.  I'm 
thinking about something like this:

start_constraining(const old_row_dofs,
                    const old_col_dofs,
                    new_row_dofs,
                    new_col_dofs);

constrain_vector(const old_row_dofs,
                  const new_row_dofs,
                  vector);

constrain_matrix(const old_row_dofs,
                  const old_col_dofs,
                  const new_row_dofs,
                  const new_col_dofs,
                  matrix);

constrain_dyad_matrix(const old_row_dofs,
                       const old_col_dofs,
                       const new_row_dofs,
                       const new_col_dofs,
                       v,
                       w);

Or, perhaps better, have a "Constraining" class that is created 
locally by the user and holds all required information (including the 
constraint matrix).

Any opinions?

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Thu, 23 Apr 2009, Tim Kroeger wrote:

> On Wed, 22 Apr 2009, Roy Stogner wrote:
>
>> One question first: can anyone think of a better name than 
>> constrain_nothing()?  The analogy to our other constrain_* functions makes 
>> sense, but it seems odd outside that context to have a method essentially 
>> named "do_nothing" that does modify its inputs.
>
> Well, of course you are right.  But actually, I find this is only a symptom 
> of the fact that a user might not expect the constrain_*() methods to modify 
> their dof_indices argument.

Good point.

> At least, that was true for me quite a long time, and it caused a
> number of programming errors in my applications.  I learned this
> when I implemented the constrain_dyad_matrix() function.  Since that
> time I understand that these methods *have* to do this.

Well, they have to get an expanded dof_indices vector.  They don't
necessarily have to expand their argument rather than making a copy or
keeping new indices in aseparate vector.  Presumably we just chose the
most efficient semantic over the more intuitive ones.

> What I want to say is this: First, having a method called
> "constrain_nothing()" might make the innocent user look into the
> details earlier and prevent him from making mistakes.

constrain_nothing() it is, then.

> Second, the whole constaining mechanism could be reworked
> completely.  I'm thinking about something like this:
>
> start_constraining(const old_row_dofs,
>                   const old_col_dofs,
>                   new_row_dofs,
>                   new_col_dofs);
>
> constrain_vector(const old_row_dofs,
>                 const new_row_dofs,
>                 vector);
>
> constrain_matrix(const old_row_dofs,
>                 const old_col_dofs,
>                 const new_row_dofs,
>                 const new_col_dofs,
>                 matrix);
>
> constrain_dyad_matrix(const old_row_dofs,
>                      const old_col_dofs,
>                      const new_row_dofs,
>                      const new_col_dofs,
>                      v,
>                      w);
>
> Or, perhaps better, have a "Constraining" class that is created locally by 
> the user and holds all required information (including the constraint 
> matrix).
>
> Any opinions?

Probably a Constraining class; otherwise we'd have to regenerate that
constraint matrix over and over again, right?  And while a more
intuitive constraint API would be nice, we've got something that works
now, so it's not a high priority for me.  Patches would be welcomed
and (eventually...) included, though.
---
Roy

Dear Roy,

On Thu, 23 Apr 2009, Roy Stogner wrote:

> constrain_nothing() it is, then.

Okay, thank you.

>> Second, the whole constaining mechanism could be reworked
>> completely.  [...]
>
> Probably a Constraining class; otherwise we'd have to regenerate that
> constraint matrix over and over again, right?

Yes, you are right.

> And while a more
> intuitive constraint API would be nice, we've got something that works
> now, so it's not a high priority for me.

The same applies to me.

> Patches would be welcomed and (eventually...) included, though.

I'll keep that in mind for the unlikely case that I should some day 
feel that I have nothing to do.  (-:

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Wed, 22 Apr 2009, Tim Kroeger wrote:

> On Fri, 17 Apr 2009, Roy Stogner wrote:
>
>> I did have time to set up your test program and get it started this 
>> morning.  Not sure if/when it's going to finish, though.
>
> Did it reveal anything yet?

On the big machine I get the same error, and relatively quickly.  I'm
going to glance over the data structures in the debugger tonight, but
unless I can find an obvious underlying cause I won't have much time
to spend on it for a while.
---
Roy

On Fri, 24 Apr 2009, Roy Stogner wrote:

> On Wed, 22 Apr 2009, Tim Kroeger wrote:
>
>> On Fri, 17 Apr 2009, Roy Stogner wrote:
>> 
>>> I did have time to set up your test program and get it started this 
>>> morning.  Not sure if/when it's going to finish, though.
>> 
>> Did it reveal anything yet?
>
> On the big machine I get the same error, and relatively quickly.  I'm
> going to glance over the data structures in the debugger tonight, but
> unless I can find an obvious underlying cause I won't have much time
> to spend on it for a while.

The old version of gdb on that machine is killing me - even with
METHOD=dbg it can't seem to find many libMesh methods, can't seem to
call the ones it does find, and sometimes crashes or returns incorrect
results from calls it does make.  If anyone knows how to walk through
the libstdc++ std::map data structure to find a specific entry, let me
know.  This didn't work for me:
http://help.lockergnome.com/linux/GDB-capabilities-exploring-STL-classes--ftopict279673.html

I've hassled the sysadmins to try and get a newer gdb to use.  For
now, examining what data structures I could directly, I've at least
got a start on the problem.

It's the same symptom as last time (constraint application failing
because a constraining DoF isn't semilocal) but it's not the same
cause.  Last time we were trying to satisfy the correct constraint
equation but we missed adding the constraining DoF to the send_list.
This time the constraint equation itself apparently came out wrong -
it's trying to use one incorrect DoF index.  And that's probably not
the root of the problem.  I presume that the odd choice of "24.01" in
the build_cube parameters was necessary to reproduce the problem?
Then the real bug has to be in one of the hacks where we use nodal
coordinates to identify a node.  I think most of my own sins in that
regard are limited to ParallelMesh.  Ben does a bit of that in
MeshCommunication that we might look at.  But the most likely culprit
is probably find_neighbors().

I've got to run right now, and won't get a chance to really look at
this again until next Wednesday.  Sorry about the delays.

And thanks again for all the ghosted vectors work.  Regardless of how
long it takes us to make it efficient, it's already been invaluable as
a debugging tool.
---
Roy

Dear Roy,

On Fri, 1 May 2009, Roy Stogner wrote:

> I presume that the odd choice of "24.01" in
> the build_cube parameters was necessary to reproduce the problem?

Actually, for me there is no reason to believe this.  The "24.01" is 
used for a different reason, having to do with some of the interna of 
the application.  I have not tested whether the crash depends on this.

(In fact, for the previous bug, as you might remember, I mixed up the 
arguments of build_cube() once again, hence unwittingly using a 
completely different geometry than in the application, which did not 
influence the bug.  Only because this time I got the argument order 
right and hence use the same geometry as in the application, this is 
not a reason to believe that the bug *depends* on the geometry.)

> I've got to run right now, and won't get a chance to really look at
> this again until next Wednesday.  Sorry about the delays.

Doesn't matter; it's not so urgent any more as it was some time ago. 
(-:

> And thanks again for all the ghosted vectors work.  Regardless of how
> long it takes us to make it efficient, it's already been invaluable as
> a debugging tool.

I am pleased to hear that!

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy,

On Mon, 4 May 2009, Tim Kroeger wrote:

>> I've got to run right now, and won't get a chance to really look at
>> this again until next Wednesday.  Sorry about the delays.
>
> Doesn't matter; it's not so urgent any more as it was some time ago.
> (-:

Well, to prevent any possible misunderstanding, I would like to add 
that by this sentence, I didn't mean that it became totally irrelevant 
for me.  Although there is currently no short term deadline pressure 
for me on this item, I would appreciate to have the ghosted vectors 
working before the next short term deadline pressure arises.  If there 
is any sensible task that I could do to assist you in finding the bug, 
please let me know.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Sat, 9 May 2009, Lorenzo Botti wrote:

> This code produces the problem during the reinit after the second solve.Hope
> it can help. 
>
> Without coarsening it seems that all works fine!

Here's what I get running the code on the libMesh svn head with METHOD=dbg:

*** Warning, This Code is Deprecated! src/base/libmesh.C, line 356, compiled May 14 2009 at 11:28:14 ***
   Beginning Solve 0
Number of elements: 219
  assembling elliptic dg system... done
System has: 768 degrees of freedom.
Linear solver converged at step: 21, final residual: 4.1012e-12
L2-Error is: 0.00666744
H1-Error is: 0.144821
   Beginning Solve 1
Number of elements: 827
  assembling elliptic dg system... done
System has: 2896 degrees of freedom.
Linear solver converged at step: 30, final residual: 1.0506e-11
L2-Error is: 0.00264921
H1-Error is: 0.102276
   Beginning Solve 2
Number of elements: 3003
  assembling elliptic dg system... done
System has: 10512 degrees of freedom.
Linear solver converged at step: 50, final residual: 2.03714e-11
L2-Error is: 0.0016323
H1-Error is: 0.070119
*** Warning, This Code is Deprecated! src/base/libmesh.C, line 366, compiled May 14 2009 at 11:28:14 ***

Since your code hadn't updated the libMesh::init() calls to use a
LibMeshInit object instead, is it safe for me to assume you're running
with an older libMesh version?  Which one?  Would you try checking out
the current SVN version and see if you can reproduce the problem
there?
---
Roy

On Thu, 14 May 2009, Tim Kroeger wrote:

> Although there is currently no short term deadline pressure for me on this 
> item, I would appreciate to have the ghosted vectors working before the next 
> short term deadline pressure arises.  If there is any sensible task that I 
> could do to assist you in finding the bug, please let me know.

If I could think of anything, I'd have mentioned it now, I'm afraid.
Right now the debug cycle is glacial, but I'm not sure how to avoid:

Recompiling.  I can't get gdb or idb to walk though our DofConstraints
structure properly (can't find function X, can't cast to unknown type
Y, etc...) and that's left me using std::cerr as a major debugging
tool.

Rerunning.  This bug depends on the precise mesh partitioning, which
we redo whenever we load a new file, so I can't just save the failing
mesh and restart from that, I have to restart from the beginning and
walk though all the dozens of AMR/C steps... in dbg or devel mode, if
I want to be able to use gdb at all on the result.

Work.  Your typical debug cycle doesn't include "Spend days giving and
listening to talks" or "Fix and rerun sensitivity analyses with a
different code", but the last few weeks have been swamped for me.

It looks like add_constraints_to_send_list didn't quite do what it was
supposed to, because this is another version of the same bug: a
hanging face node has four dependencies, one of which is a hanging
edge node with two dependencies, and the processor with the face node
somehow isn't getting the farther grand-dependency added to its
send_list.

It'll be another run or two until I'm certain of that, though, and
probably a few more before I've figured out why it's happening and
how to fix it.
---
Roy

Dear Roy,

On Thu, 14 May 2009, Roy Stogner wrote:

> Right now the debug cycle is glacial, but I'm not sure how to avoid: 
> [...]

Thank you very much for your interim report.  I just wanted to make 
sure that things have not been forgotten.  I understand that it takes 
time to find this bug.

Sorry if my mail made you feel defensive.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Fri, 15 May 2009, Tim Kroeger wrote:

> On Thu, 14 May 2009, Roy Stogner wrote:
>
>> Right now the debug cycle is glacial, but I'm not sure how to avoid: [...]
>
> Thank you very much for your interim report.  I just wanted to make sure that 
> things have not been forgotten.  I understand that it takes time to find this 
> bug.
>
> Sorry if my mail made you feel defensive.

Not defensive at you, just exasperated at gdb.  My attempts to get
a proper debugger's-eye-view of complex STL-tree-based classes have
been both numerous and fruitless.
---
Roy

Roy Stogner wrote:

> My attempts to get a proper debugger's-eye-view of complex
> STL-tree-based classes have been both numerous and fruitless.

Do you have any experience with

  http://sourceware.org/gdb/wiki/ProjectArcher

I haven't used it because I work primarily in C, but it should help with
this.

Jed

On Fri, 15 May 2009, Jed Brown wrote:

> Roy Stogner wrote:
>
>> My attempts to get a proper debugger's-eye-view of complex
>> STL-tree-based classes have been both numerous and fruitless.
>
> Do you have any experience with
>
>  http://sourceware.org/gdb/wiki/ProjectArcher

No, thank you!  I'd tried a couple different macro sets that were
supposed to work on top of vanilla gdb, but never anything that
changed the source itself.  I'll give this a shot.
---
Roy

On Thu, 14 May 2009, Roy Stogner wrote:

> It'll be another run or two until I'm certain of that, though, and
> probably a few more before I've figured out why it's happening and
> how to fix it.

Got it!

The reason this bug was so hard to find is that, in some sense, it's
more of a design mistake than a bug!  enforce_constraints_exactly()
works correctly with serial vectors; ghosted vectors are correctly
behaving as their API specifies... the two correct behaviors just
weren't 100% compatible yet, is all.

When I wrote enforce_constraints_exactly(), for some reason, (possibly
because we weren't yet recursively expanding contraints and so a more
straightforward method wasn't yet possible, possibly because I was
still learning how the contraint system worked), I made it work by
building the constraint matrix C for each element, looping over the
rows that correspond to local constrained dofs, and setting 
vglobal_i = sum(C_ij*vlocal_j)

If you'll forgive my abused notation: The problem here is that we've
got an element with dof a on processor 1 and dof b on processor 2, the
latter of which depends on dof c on processor 3.  Because dof c isn't
a constraint dependency on processor 1, that processor doesn't have it
in the send_list.  This means that vlocal_c is inaccurately 0 on a
processor 1 serial vector, but we don't really care, because C_ac is 0
too and the inaccuracy in vlocal_c doesn't propagate to vglobal_a.
But accessing vlocal_c, even to multiply it by 0, throws a
libmesh_error() on a ghosted vector!

Anyway, an immediate fix is simple: just skip accumulating indices
where C_ij==0.0.  I've committed that to SVN, and on my machine it
works to take the test case you sent all the way to completion.  In
the long term we'll want to change enforce_constraints_exactly to just
loop over local dofs and directly use the constraint rows, but I'd
like to make sure this old bug is fixed before I risk mucking things
up and adding new bugs.  ;-)
---
Roy

On Sat, May 16, 2009 at 8:26 AM, Roy Stogner <roystgnr@...> wrote:
>
> On Thu, 14 May 2009, Roy Stogner wrote:
>
>> It'll be another run or two until I'm certain of that, though, and
>> probably a few more before I've figured out why it's happening and
>> how to fix it.
>
> Got it!
>
> The reason this bug was so hard to find is that, in some sense, it's
> more of a design mistake than a bug!  enforce_constraints_exactly()
> works correctly with serial vectors; ghosted vectors are correctly
> behaving as their API specifies... the two correct behaviors just
> weren't 100% compatible yet, is all.
>
> When I wrote enforce_constraints_exactly(), for some reason, (possibly
> because we weren't yet recursively expanding contraints and so a more
> straightforward method wasn't yet possible, possibly because I was
> still learning how the contraint system worked), I made it work by
> building the constraint matrix C for each element, looping over the
> rows that correspond to local constrained dofs, and setting
> vglobal_i = sum(C_ij*vlocal_j)
>
> If you'll forgive my abused notation: The problem here is that we've
> got an element with dof a on processor 1 and dof b on processor 2, the
> latter of which depends on dof c on processor 3.  Because dof c isn't
> a constraint dependency on processor 1, that processor doesn't have it
> in the send_list.  This means that vlocal_c is inaccurately 0 on a
> processor 1 serial vector, but we don't really care, because C_ac is 0
> too and the inaccuracy in vlocal_c doesn't propagate to vglobal_a.
> But accessing vlocal_c, even to multiply it by 0, throws a
> libmesh_error() on a ghosted vector!
>
> Anyway, an immediate fix is simple: just skip accumulating indices
> where C_ij==0.0.  I've committed that to SVN, and on my machine it
> works to take the test case you sent all the way to completion.  In
> the long term we'll want to change enforce_constraints_exactly to just
> loop over local dofs and directly use the constraint rows, but I'd
> like to make sure this old bug is fixed before I risk mucking things
> up and adding new bugs.  ;-)

Way to go Roy!

I'd like to personally double your libmesh developer salary this month ;-)

-- 
John

> Way to go Roy!
> 
> I'd like to personally double your libmesh developer salary this month ;-)

I second that.  

When we are all in Austin week after next let's decide on a libMesh-0.7.0
release date.  We also need to devise some bonus incentive program for
meeting the milestone? ;-)

-Ben

Dear Roy,

On Sat, 16 May 2009, Roy Stogner wrote:

> Got it! [...]

Great work!  I'll let the application run today and see whether any 
new bugs emerge.  (-:

> Anyway, an immediate fix is simple: just skip accumulating indices
> where C_ij==0.0.  I've committed that to SVN, and on my machine it
> works to take the test case you sent all the way to completion.

I would strongly suggest to add a comment at that position in the 
code.  Otherwise, if you (or somebody else) later should decide that 
-Wfloat-equal should be enabled, you'll get a warning at this point, 
and since you might have forgotten the reason, you'll be tempted to 
remove that seemingly useless piece of code.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy,

On Mon, 18 May 2009, Tim Kroeger wrote:

>> Got it! [...]
>
> Great work!  I'll let the application run today and see whether any
> new bugs emerge.  (-:

When I first tested it, it kept crashing at the same point as before. 
It took me quite long to find out what I did wrong: I compiled using

 	METHO=devel make

and linked against the devel version.  Seems right, doesn't it?  Well, 
but if you look closely, you see that I missed out a "D", so I 
actually compiled the optimized version.

Feature request for your build system: If the user ever has compiled a 
version other than optimized, "make" with unset "METHOD" variable 
should produce an error message and not work.

After I had fixed my compile statement, the next run crashed at a 
different (later) point.  I conjecture that that was out-of-memory, so 
I restarted it on a larger number of nodes (same number of CPUs, 
though).  It then finally ran through without crash.

By the way, the temporal scalability is also quite bad.  You might 
remember that this bug has emerged as a test of how my application 
behaves on a larger number of CPUs.  The runtimes for 8 CPUs on 3 
nodes (copied from my mail of April 8, 2009) were:

no-ghosted-1 : 11:34:10
no-ghosted-2 : 11:35:54
ghosted-1    : 17:25:28
ghosted-2    : 16:33:23

Now, the new result for 24 CPUs on 4 nodes (3 nodes ran out of memory) 
is:

ghosted      : 10:42:34

Not really nice.  In particular not essentially faster than the 
unghosted version on a comparable number of nodes.

I'll perfom the log output that Jed suggested and see what he says. 
Perhaps there is some easy possibility to make it faster.

>> Anyway, an immediate fix is simple: just skip accumulating indices
>> where C_ij==0.0.  I've committed that to SVN, and on my machine it
>> works to take the test case you sent all the way to completion.
>
> I would strongly suggest to add a comment at that position in the
> code.  Otherwise, if you (or somebody else) later should decide that
> -Wfloat-equal should be enabled, you'll get a warning at this point,
> and since you might have forgotten the reason, you'll be tempted to
> remove that seemingly useless piece of code.

I guess that you didn't have the time yet to write such a comment, so 
I wrote that for you (see attachment), you just have to check it in.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:

> Not really nice.  In particular not essentially faster than the
> unghosted version on a comparable number of nodes.

It is important to know what preconditioners are being used
(preconditioners always change in parallel, though not as much when
there is a coarse level as in multigrid).  Also, memory performance
(especially bandwidth) is usually the overwhelming issue for implicit
solvers (e.g. you are very lucky to get 4% of peak FPU performance for
MatVec on Core 2 Quad).  Thus using more cores frequently does not help,
and you need more sockets to improve performance.  Network latency is
also a factor, but if your subdomains are big enough that you are having
memory issues, it should not be an issue (rather, if it is an issue then
we have faulty algorithms at play).

There may well be more going on here, but it's important to consider
these issues.

> I'll perfom the log output that Jed suggested and see what he says.
> Perhaps there is some easy possibility to make it faster.

FWIW, I always run nontrivial jobs with -log_summary.  It does not add
measurable overhead and that output is really useful.  If the effort
required to ignore that output is less than the run (i.e. the run is
more than a few seconds), it is worthwhile to have that profiling info.

Jed

Dear Jed,

On Tue, 26 May 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>
>> Not really nice.  In particular not essentially faster than the
>> unghosted version on a comparable number of nodes.
>
> It is important to know what preconditioners are being used
> (preconditioners always change in parallel, though not as much when
> there is a coarse level as in multigrid).

Seems as if I'm using BLOCK_JACOBI in parallel.

Also, the preconditioner is not applied to the original matrix (it's 
the issue with the ShellMatrix, which you certainly remember), but 
this can't be the problem here because the matrix is independent of 
the number of CPUs.

> Also, memory performance
> (especially bandwidth) is usually the overwhelming issue for implicit
> solvers (e.g. you are very lucky to get 4% of peak FPU performance for
> MatVec on Core 2 Quad).  Thus using more cores frequently does not help,
> and you need more sockets to improve performance.

I must admit that I'm not familiar with these hardware issues at all, 
so your explanations are all Greek to me.  Well, since you were 
talking about sockets, I tried this (on one of the nodes):

tkroeger@...:~> cat /proc/net/sockstat
sockets: used 239
TCP: inuse 31 orphan 0 tw 21 alloc 39 mem 2
UDP: inuse 19
RAW: inuse 0
FRAG: inuse 0 memory 0

Is this the kind of information that can help finding what is going 
on?  If not, please let me know how I can find the required 
information.

>> I'll perfom the log output that Jed suggested and see what he says.
>> Perhaps there is some easy possibility to make it faster.
>
> FWIW, I always run nontrivial jobs with -log_summary.  It does not add
> measurable overhead and that output is really useful.

I've enabled it now and started the application.  I have not yet 
reduced the complexity (number of time steps) because, for the moment, 
it doesn't bother me if I get the result on Monday.  If it should turn 
out that more tests have to be made, I will think about that again.

(BTW, in the last two weeks, I was very busy with some other work that 
had a tight deadline, but now I should have more time again to work on 
this issue.)

Anyway, I'll send you the output of -log_summary as soon as I have it.

Thank you very much in advance for your assistance!

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:
> Dear Jed,
> 
> On Tue, 26 May 2009, Jed Brown wrote:
> 
>> Tim Kroeger wrote:
>>
>>> Not really nice.  In particular not essentially faster than the
>>> unghosted version on a comparable number of nodes.
>>
>> It is important to know what preconditioners are being used
>> (preconditioners always change in parallel, though not as much when
>> there is a coarse level as in multigrid).
> 
> Seems as if I'm using BLOCK_JACOBI in parallel.

This is the default and it changes a lot when you change the number of
processors.  You could easily be needing a lot more iterations in the
linear solve when you use more processors.  Comparing -log_summary
output would tell us if this is so.  Unless you use a direct solve or
have multilivel domain decomposition (e.g. multigrid or iterative
substructuring like FETI-DP or BDDC) working for you, the number of
iterations will increase when you add processors.  There is an
asymptotic theory, but in practice this effect is highly
problem-dependent.

> Also, the preconditioner is not applied to the original matrix (it's the
> issue with the ShellMatrix, which you certainly remember), but this
> can't be the problem here because the matrix is independent of the
> number of CPUs.

This is normal, I never have the same Jacobian as preconditioning
matrix.

>> Also, memory performance
>> (especially bandwidth) is usually the overwhelming issue for implicit
>> solvers (e.g. you are very lucky to get 4% of peak FPU performance for
>> MatVec on Core 2 Quad).  Thus using more cores frequently does not help,
>> and you need more sockets to improve performance.
> 
> I must admit that I'm not familiar with these hardware issues at all, so
> your explanations are all Greek to me.  Well, since you were talking
> about sockets, I tried this (on one of the nodes):
> 
> tkroeger@...:~> cat /proc/net/sockstat
> sockets: used 239
> TCP: inuse 31 orphan 0 tw 21 alloc 39 mem 2
> UDP: inuse 19
> RAW: inuse 0
> FRAG: inuse 0 memory 0
> 
> Is this the kind of information that can help finding what is going on? 
> If not, please let me know how I can find the required information.

Heh, I was referring to hardware sockets, not network sockets.

So most CPUs these days are multicore, but the cores share the same
memory bus (at least partially).  When you say 8 CPUs on 3 nodes or 24
CPUs on 4 nodes, I assume you mean 8 cores on 3 nodes, 24 cores on 4
nodes.  This is a different number of cores per node, but the
configuration of these nodes is relevant.  For instance, 4 dual-core is
almost always faster than 2 quad-core (although the latter is probably
cheaper and uses less power).  Sparse matrix kernels are typically much
faster per process if there is only one process per socket (i.e. per
CPU, not per core).  The floating point unit (FPU) on a single core is
so much faster than memory that even if you only use one core (on Core
2) you will be happy to see 10% of peak FPU throughput.  If you use all
4 cores and have custom sparse matrix kernels, you might see 4% of peak.
Now 4*4% = 16% is more than 10% so you are getting more work done by
using those extra cores, but nowhere near four times more work.

Getting 4% requires a huge amount of effort.  In particular, threads
need to be aware of what the other cores have in L1 cache.  This doesn't
work in MPI and is not portable, so you'll actually see somewhat less
performance when using any portable library.  This is why CPU speed is
nearly irrelevant for sparse linear algebra and you are much better off
to buy memory bandwidth.  (Explicit methods, residual evaluation, and
Jacobian assembly often benefit from a faster CPU so it's not always a
waste to buy a fast CPU, it just doesn't matter much inside the linear
solver.)

Knowing your hardware (/proc/cpuinfo and network) and MPI configuration
(which MPI and which device it is using) is useful.

Jed

Dear Jed,

On Fri, 5 Jun 2009, Jed Brown wrote:

>>> Also, memory performance
>>> (especially bandwidth) is usually the overwhelming issue for implicit
>>> solvers (e.g. you are very lucky to get 4% of peak FPU performance for
>>> MatVec on Core 2 Quad).  Thus using more cores frequently does not help,
>>> and you need more sockets to improve performance.
>>
>> I must admit that I'm not familiar with these hardware issues at all, so
>> your explanations are all Greek to me.  Well, since you were talking
>> about sockets, I tried this (on one of the nodes):
>>
>> tkroeger@...:~> cat /proc/net/sockstat
>> sockets: used 239
>> TCP: inuse 31 orphan 0 tw 21 alloc 39 mem 2
>> UDP: inuse 19
>> RAW: inuse 0
>> FRAG: inuse 0 memory 0
>>
> Heh, I was referring to hardware sockets, not network sockets.

Well, I told you I have no idea what I'm talking about.  (-:

> So most CPUs these days are multicore, but the cores share the same
> memory bus (at least partially).  When you say 8 CPUs on 3 nodes or 24
> CPUs on 4 nodes, I assume you mean 8 cores on 3 nodes, 24 cores on 4
> nodes.  This is a different number of cores per node, but the
> configuration of these nodes is relevant.  For instance, 4 dual-core is
> almost always faster than 2 quad-core (although the latter is probably
> cheaper and uses less power).

/proc/cpuinfo is attached.  Seems as if I've got 2 quad-core CPUs on 
each node, right?

> Knowing your hardware (/proc/cpuinfo and network) and MPI configuration
> (which MPI and which device it is using) is useful.

Hmm, how do I find out about the network hardware?

There is a pre-installed MPI version, that is:

ParaStation 5.0.9-1 (Wed Nov 12 18:27:53 CET 2008)
  config args: '--prefix=/opt/parastation/mpi2' '--with-confset=default'
  mpich2 version: 1.0.7

But for some cumbersome reasons I can't use it.  Well, perhaps I just 
haven't tried hard enough.  I installed another MPI myself, and that 
is mpich2 1.0.8, without any options (other than the installation 
path).  I'll send you its config.log off-list (180KB).  Well, now that 
I think about it, I wonder whether I could simply supply the above 
stated config arguments manually to my local MPI installation.  Would 
that be a promising idea?

Also, it might be useful to note that the cluster has two different 
queues, one of which is called "ethernet" and the other "infiniband". 
They use different nodes.  These names seem to indicate some 
information about the hardware (let's assume they have been chosen on 
purpose).  I understand that "infiniband" refers to the "better" 
hardware in some sense, but the "infiniband" queue has a much smaller 
number of nodes.  Also, the admin told me that using my manually 
installed MPI, I would not be able to benefit from the advantages that 
the "inifiniband" queue has.  Hence, I was using the "ethernet" queue 
all the time.  /proc/cpuinfo does not differ between the nodes of the 
two queues (except for changes in some minor decimal places of the 
values of "cpu MHz" and "bogomips").

Anyway, I'll send you my log output on Monday.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:
> There is a pre-installed MPI version, that is:
> 
> ParaStation 5.0.9-1 (Wed Nov 12 18:27:53 CET 2008)
>  config args: '--prefix=/opt/parastation/mpi2' '--with-confset=default'
>  mpich2 version: 1.0.7
> 
> But for some cumbersome reasons I can't use it.  Well, perhaps I just
> haven't tried hard enough.  I installed another MPI myself, and that is
> mpich2 1.0.8, without any options (other than the installation path). 
> I'll send you its config.log off-list (180KB).  Well, now that I think
> about it, I wonder whether I could simply supply the above stated config
> arguments manually to my local MPI installation.  Would that be a
> promising idea?

I don't know what confset is, but your install is probably fine.  So the
default device through 1.0.8 was ch3:sock which uses sockets (the
low-level API for moving information over the network, these are the
sockets that sockstat referred to ;-) for all communication.  This is
fine, but you would probably get better performance with ch3:nemesis
(default now in 1.1, released on Wed) which uses shared memory locally
and sockets (or Myrinet if you have it) to talk to remote processes.

> Also, it might be useful to note that the cluster has two different
> queues, one of which is called "ethernet" and the other "infiniband".
> They use different nodes.  These names seem to indicate some information
> about the hardware (let's assume they have been chosen on purpose).  I
> understand that "infiniband" refers to the "better" hardware in some
> sense, but the "infiniband" queue has a much smaller number of nodes. 
> Also, the admin told me that using my manually installed MPI, I would
> not be able to benefit from the advantages that the "inifiniband" queue
> has.  Hence, I was using the "ethernet" queue all the time. 
> /proc/cpuinfo does not differ between the nodes of the two queues
> (except for changes in some minor decimal places of the values of "cpu
> MHz" and "bogomips").

MPICH2 doesn't support IB, you would need MVAPICH2 or Open MPI or a
vendor-supplied MPI (likely to be a derivative of one of these).  With
gigabit Ethernet, we expect latencies around 100 microseconds which is a
lot more than IB (< 2us).  Depending on the decomposition and size of
the subdomains, this latency could already be an issue.  But since you
were having issues running out of memory, I'm assuming that your
subdomains are pretty big in which case network latency shouldn't be
crippling (there is an algorithmic bug if it is).

Do you have any of the usual libmesh profiling output?  The granularity
is pretty coarse, but it would indicate if most of your time is spent in
the linear solve or assembly or whatever.

Jed

Dear Jed,

To start with the bad news: There is no log output yet.  The 
application hasn't finished.  It reduced the time step size to 2e-6 
for some reason.  I can't understand this because I'm not aware of any 
changes I made since the last run -- except for enabling 
"-log_summary", but that should not change the behaviour, should it?

I will try to find out today what happend.

On Fri, 5 Jun 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>> There is a pre-installed MPI version, that is:
>>
>> ParaStation 5.0.9-1 (Wed Nov 12 18:27:53 CET 2008)
>>  config args: '--prefix=/opt/parastation/mpi2' '--with-confset=default'
>>  mpich2 version: 1.0.7
>>
>> But for some cumbersome reasons I can't use it.  Well, perhaps I just
>> haven't tried hard enough.  I installed another MPI myself, and that is
>> mpich2 1.0.8, without any options (other than the installation path).
>> I'll send you its config.log off-list (180KB).  Well, now that I think
>> about it, I wonder whether I could simply supply the above stated config
>> arguments manually to my local MPI installation.  Would that be a
>> promising idea?
>
> I don't know what confset is, but your install is probably fine.

Good to hear.

> So the
> default device through 1.0.8 was ch3:sock which uses sockets (the
> low-level API for moving information over the network, these are the
> sockets that sockstat referred to ;-) for all communication.  This is
> fine, but you would probably get better performance with ch3:nemesis
> (default now in 1.1, released on Wed) which uses shared memory locally
> and sockets (or Myrinet if you have it) to talk to remote processes.

So you would recommend to update to mpich2 1.1?

(If I do, I think I will have to recompile PETSc, libmesh, and my 
application as well.)

>> Also, it might be useful to note that the cluster has two different
>> queues, one of which is called "ethernet" and the other 
>> "infiniband". [...]
>
> MPICH2 doesn't support IB, you would need MVAPICH2 or Open MPI or a
> vendor-supplied MPI (likely to be a derivative of one of these).  With
> gigabit Ethernet, we expect latencies around 100 microseconds which is a
> lot more than IB (< 2us).  Depending on the decomposition and size of
> the subdomains, this latency could already be an issue.  But since you
> were having issues running out of memory, I'm assuming that your
> subdomains are pretty big in which case network latency shouldn't be
> crippling (there is an algorithmic bug if it is).

Okay, so I will continue using the "ethernet" queue.

> Do you have any of the usual libmesh profiling output?

No, I haven't at the moment.  Does it suffice to do "--enable-perflog" 
as a configure option for that?  I'll give it a try.

But first, I think I should find out why the application behaves so 
much different now from what it did a couple of weeks ago.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Jed,

On Mon, 8 Jun 2009, Tim Kroeger wrote:

> To start with the bad news: There is no log output yet.  The
> application hasn't finished.  It reduced the time step size to 2e-6
> for some reason.  I can't understand this because I'm not aware of any
> changes I made since the last run -- except for enabling
> "-log_summary", but that should not change the behaviour, should it?

Today's good news: It's running "correctly" now.

Today's fair news: No log output yet, but to be expected tomorrow.

Today's bad news: The term "correctly" in the above sentence means 
"wrong".  The problem that I stated yesterday was due to a bugfix I 
performed in the meantime.  I reverted that fix now, but without that 
fix, the result that my application gets is likely to not have 
anything to do with the thing it was supposed to compute.  Actually, I 
can't understand why the application keeps producing results that look 
sensible.

Anyway, I decided to postpone that problem for a while and concentrate 
now on the performance issues.

I haven't updated to mpich2 1.1 yet; let's leave that as an option for 
later.  However, I did --enable-perflog and -log_summary.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Jed,

On Tue, 9 Jun 2009, Tim Kroeger wrote:

> Today's fair news: No log output yet, but to be expected tomorrow.

Please find attached the log output now.

As I see it, there are two obvious things that I could do.  That is 
(1) update to mpich2 1.1, and (2) recompile PETSc with 
"--with-debugging=no".  Would you recommend to do so?

On the other hand, a main issue was that the ghosted code is by about 
a factor of 1.5 slower than the non-ghosted code.  That shouldn't have 
to do with debugging options, because they were enabled for both 
versions.

Let me know what you think and what you can read from the logs.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Jed and Lorenzo,

On Wed, 10 Jun 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>
>> Please find attached the log output now.
>
> First, I'm confused about why all the detailed output is missing.  See
> attached for the sort of output I expect.  In particular, there is
> fine-grained output on vector operations and solver internals with
> details about load balancing.

I have no idea what happend here.  I should note that I supplied the 
option "-log_summary" using PetscOptionsSetValue rather than directly 
on the command line, but that shouldn't be the point, should it?

> I see about 30% of the time spent in assemble which about what I'd
> expect, but solve time is relatively low, especially when compared to
> some of the FE ops.  I don't recall, is a significant part of your code
> explicit?

There are some explicit parts, but it wouldn't expect them to be 
significant.

When I watch the output while the application is running, it feels as 
if it is "most of the time" either inside assemble/solve of one of the 
implicit parts or inside EquationSystems::reinit().

> inverse_map is getting called a ton.

Yes.  I use a PointLocator, but that is at then end, after the main 
computation has finished.  The PointLocator is used very heavily then, 
and this part of the code runs in serial (only on core #0) and takes 
about 35 minutes.

I see that this contradicts the fact that over 90 minutes are spent 
inside inverse_map().  I have currently no idea where else I am using 
that.

>> As I see it, there are two obvious things that I could do.  That is (1)
>> update to mpich2 1.1, and (2) recompile PETSc with
>> "--with-debugging=no".  Would you recommend to do so?
>
> 1. Time for MPI_Barrier/MPI_Send looks a bit high, that might be
> improved with ch3:nemesis, but would mainly be due to network
> configuration.
>
> 2. This can make a significant difference on the PETSc operations, but
> it's not clear that much time is being spent in PETSc anyway.  If you do
> rebuild, maybe you could also upgrade to PETSc-3.0.0.

Okay, so I won't rebuild now but keep this in mind.

>> On the other hand, a main issue was that the ghosted code is by about a
>> factor of 1.5 slower than the non-ghosted code.  That shouldn't have to
>> do with debugging options, because they were enabled for both versions.
>
> Do you have a log from an equivalent non-ghosted run?

Good point.  (-:

I don't have one right now.  Actually, running on 24 cores on 4 nodes 
won't be possible without ghosted (due to memory shortage).  I will 
perform two computations (one ghosted and one non-ghosted) with 8 
cores on 4 nodes and send you the logs.  That will probably be 
tomorrow.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:
> Dear Jed and Lorenzo,
> 
> On Wed, 10 Jun 2009, Jed Brown wrote:
> 
>> Tim Kroeger wrote:
>>
>>> Please find attached the log output now.
>>
>> First, I'm confused about why all the detailed output is missing.  See
>> attached for the sort of output I expect.  In particular, there is
>> fine-grained output on vector operations and solver internals with
>> details about load balancing.
> 
> I have no idea what happend here.  I should note that I supplied the
> option "-log_summary" using PetscOptionsSetValue rather than directly on
> the command line, but that shouldn't be the point, should it?

This is not okay, that just registers the option in the options
database, but the profiling functions don't check there (for performance
reasons).  The option either needs to be present at PetscInitialize,
i.e. any of

1. Passed in on the command-line
2. Specified in the file that is the third argument of PetscInitialize
3. Put in $HOME/.petscrc
4. Put in $PWD/.petscrc

OR, if you must, call PetscLogBegin() to start logging and make sure
that -log_summary is in the options database by the time PetscFinalize
is called.

>> I see about 30% of the time spent in assemble which about what I'd
>> expect, but solve time is relatively low, especially when compared to
>> some of the FE ops.  I don't recall, is a significant part of your code
>> explicit?
> 
> There are some explicit parts, but it wouldn't expect them to be
> significant.
> 
> When I watch the output while the application is running, it feels as if
> it is "most of the time" either inside assemble/solve of one of the
> implicit parts or inside EquationSystems::reinit().

According to the libmesh profiling, the sum of these steps is only half
the runtime.  It should be essentially everything (especially since
function evaluation is not a separate step).  But, it appears that the
final profiling output for, say, assemble(), does not include sub-events
(so there is no way for us to know the total amount of time spent in
"assembly" without defining PerfLog objects for the various stages).
This is different than PETSc's convention so keep that in mind.

>> inverse_map is getting called a ton.
> 
> Yes.  I use a PointLocator, but that is at then end, after the main
> computation has finished.  The PointLocator is used very heavily then,
> and this part of the code runs in serial (only on core #0) and takes
> about 35 minutes.
> 
> I see that this contradicts the fact that over 90 minutes are spent
> inside inverse_map().  I have currently no idea where else I am using that.

It's also called in assembly somewhere and seems to be a significant
bottleneck.  Roy et al., could you remind me why this is required?

>>> As I see it, there are two obvious things that I could do.  That is (1)
>>> update to mpich2 1.1, and (2) recompile PETSc with
>>> "--with-debugging=no".  Would you recommend to do so?
>>
>> 1. Time for MPI_Barrier/MPI_Send looks a bit high, that might be
>> improved with ch3:nemesis, but would mainly be due to network
>> configuration.
>>
>> 2. This can make a significant difference on the PETSc operations, but
>> it's not clear that much time is being spent in PETSc anyway.  If you do
>> rebuild, maybe you could also upgrade to PETSc-3.0.0.
> 
> Okay, so I won't rebuild now but keep this in mind.
> 
>>> On the other hand, a main issue was that the ghosted code is by about a
>>> factor of 1.5 slower than the non-ghosted code.  That shouldn't have to
>>> do with debugging options, because they were enabled for both versions.
>>
>> Do you have a log from an equivalent non-ghosted run?
> 
> Good point.  (-:
> 
> I don't have one right now.  Actually, running on 24 cores on 4 nodes
> won't be possible without ghosted (due to memory shortage).  I will
> perform two computations (one ghosted and one non-ghosted) with 8 cores
> on 4 nodes and send you the logs.  That will probably be tomorrow.

Great.


Jed

On Thu, 11 Jun 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>
>> I see that this contradicts the fact that over 90 minutes are spent
>> inside inverse_map().  I have currently no idea where else I am using that.
>
> It's also called in assembly somewhere and seems to be a significant
> bottleneck.  Roy et al., could you remind me why this is required?

It's used in assembly of DG systems - when taking jump terms on
element sides we take an inverse map of the physical coordinates found
on one element to get to the master coordinates on the neighbor
element.  Likewise for CG systems using jump error estimators.

It's used in elem->contains_point() which a PointLocator or
MeshFunction would be hitting a lot.

We're using it unnecessarily and dangerously in find_point_neighbors;
that won't be an issue for Tim but I ought to fix that before I start
heavily using patch recovery estimators again...

And it's used a lot in System::project_vector, which is probably where
Tim's code is hitting it.  I'm still stunned that it would be taking
up 90 minutes, though!  On his affine hex meshes that call ought to be
converging in one step after a couple dozen flops...
---
Roy

Dear Jed,

On Thu, 11 Jun 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>>
>> I have no idea what happend here.  I should note that I supplied the
>> option "-log_summary" using PetscOptionsSetValue rather than directly on
>> the command line, but that shouldn't be the point, should it?
>
> This is not okay, that just registers the option in the options
> database, but the profiling functions don't check there (for performance
> reasons).  The option either needs to be present at PetscInitialize,
> i.e. any of
>
> 1. Passed in on the command-line
> 2. Specified in the file that is the third argument of PetscInitialize
> 3. Put in $HOME/.petscrc
> 4. Put in $PWD/.petscrc
>
> OR, if you must, call PetscLogBegin() to start logging and make sure
> that -log_summary is in the options database by the time PetscFinalize
> is called.

Ah, I didn't know that.  I will call PetscLogBegin() directly and 
restart my application.  This will postpone the logfiles until Monday 
(except for a small chance that I might go to work on the weekend).

>> When I watch the output while the application is running, it feels as if
>> it is "most of the time" either inside assemble/solve of one of the
>> implicit parts or inside EquationSystems::reinit().
>
> According to the libmesh profiling, the sum of these steps is only half
> the runtime.  It should be essentially everything (especially since
> function evaluation is not a separate step).  But, it appears that the
> final profiling output for, say, assemble(), does not include sub-events
> (so there is no way for us to know the total amount of time spent in
> "assembly" without defining PerfLog objects for the various stages).
> This is different than PETSc's convention so keep that in mind.

Yes, that's the point.  I think this was different in earlier times, 
and they changed it then, because in some cases it's more useful the 
way it is done now.  Unfortunately not for all.

I've now manipulated the PerfLog class such that it does *both*.

Roy: I have attached a patch concerning this point.  What do you 
think about it.  I would appreciate if it would be committed. 
You see, it doesn't remove functionality; however, I understand 
that it makes the output of PerfLog wider by 24 charaters and 
thus eventually less well to read.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Sorry, forgot the patch.  Here it is.

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy,

On Thu, 11 Jun 2009, Roy Stogner wrote:

> On Thu, 11 Jun 2009, Jed Brown wrote:
>
>> Tim Kroeger wrote:
>> 
>>> I see that this contradicts the fact that over 90 minutes are spent
>>> inside inverse_map().  I have currently no idea where else I am using 
>>> that.
>> 
>> It's also called in assembly somewhere and seems to be a significant
>> bottleneck.  Roy et al., could you remind me why this is required?
>
> It's used in assembly of DG systems - when taking jump terms on
> element sides we take an inverse map of the physical coordinates found
> on one element to get to the master coordinates on the neighbor
> element.  Likewise for CG systems using jump error estimators.

Okay, I'm not using those systems.  (Actually, I don't even know what 
they do.  Sometimes I think I should use them...)

> It's used in elem->contains_point() which a PointLocator or
> MeshFunction would be hitting a lot.

I do that after the computation (as described before), but that's only 
35 minutes total time.

> We're using it unnecessarily and dangerously in find_point_neighbors;
> that won't be an issue for Tim but I ought to fix that before I start
> heavily using patch recovery estimators again...

I can't say anything about this since it seems to be some libMesh 
internal thing.  Well, Roy seems to know that I don't use that, so 
it's fine.  (-:

> And it's used a lot in System::project_vector, which is probably where
> Tim's code is hitting it.  I'm still stunned that it would be taking
> up 90 minutes, though!  On his affine hex meshes that call ought to be
> converging in one step after a couple dozen flops...

Hmm.

I've now removed that part of the code in my application that uses the 
MeshFunction after the main computation (see above).  The new log 
outputs should give us more insight.  In particular, we will see 
exactly how much time is spent inside those inverse_map() calls that 
are performed in System::project_vector().

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:
> Dear Roy,
> 
> On Thu, 11 Jun 2009, Roy Stogner wrote:
> 
>> On Thu, 11 Jun 2009, Jed Brown wrote:
>>
>>> Tim Kroeger wrote:
>>>
>>>> I see that this contradicts the fact that over 90 minutes are spent
>>>> inside inverse_map().  I have currently no idea where else I am
>>>> using that.
>>>
>>> It's also called in assembly somewhere and seems to be a significant
>>> bottleneck.  Roy et al., could you remind me why this is required?
>>
>> It's used in assembly of DG systems - when taking jump terms on
>> element sides we take an inverse map of the physical coordinates found
>> on one element to get to the master coordinates on the neighbor
>> element.  Likewise for CG systems using jump error estimators.
> 
> Okay, I'm not using those systems.  (Actually, I don't even know what
> they do.  Sometimes I think I should use them...)
> 
>> It's used in elem->contains_point() which a PointLocator or
>> MeshFunction would be hitting a lot.
> 
> I do that after the computation (as described before), but that's only
> 35 minutes total time.
> 
>> We're using it unnecessarily and dangerously in find_point_neighbors;
>> that won't be an issue for Tim but I ought to fix that before I start
>> heavily using patch recovery estimators again...
> 
> I can't say anything about this since it seems to be some libMesh
> internal thing.  Well, Roy seems to know that I don't use that, so it's
> fine.  (-:
> 
>> And it's used a lot in System::project_vector, which is probably where
>> Tim's code is hitting it.  I'm still stunned that it would be taking
>> up 90 minutes, though!  On his affine hex meshes that call ought to be
>> converging in one step after a couple dozen flops...
> 
> Hmm.

It looks like it is called for every quadrature point, probably the most
expensive part of FE::reinit.

(gdb) bt
#0  FE<3u, (libMeshEnums::FEFamily)0>::inverse_map (elem=0x20f00a0, physical_point=@0x255b820, tolerance=9.9999999999999995e-07, secure=true) at src/fe/fe_map.C:635
#1  0x00007f48c7ac893e in FE<3u, (libMeshEnums::FEFamily)0>::inverse_map (elem=0x20f00a0, physical_points=@0x252dbf0, reference_points=@0x7fffd0b06bf0, 
    tolerance=9.9999999999999995e-07, secure=true) at src/fe/fe_map.C:985
#2  0x00007f48c79e9fcf in FE<3u, (libMeshEnums::FEFamily)0>::reinit (this=0x252dbe0, elem=0x20f00a0, s=0, tolerance=9.9999999999999995e-07) at src/fe/fe_boundary.C:135
#3  0x000000000044717f in assemble_poisson (es=@0x7fffd0b088f0, system_name=@0x215d578) at ex4.C:521
#4  0x00007f48c8229d6f in System::user_assembly (this=0x215d520) at src/solvers/system.C:1054
#5  0x00007f48c8223661 in System::assemble (this=0x215d520) at src/solvers/system.C:339
#6  0x00007f48c820b824 in ImplicitSystem::assemble (this=0x215d520) at src/solvers/implicit_system.C:163
#7  0x00007f48c821072b in LinearImplicitSystem::solve (this=0x215d520) at src/solvers/linear_implicit_system.C:82
#8  0x000000000044560a in main (argc=5, argv=0x7fffd0b08cf8) at ex4.C:239


But this is not dependent on whether ghosted vectors are used, so I'm
really curious to see that output.

Jed

On Thu, 11 Jun 2009, Jed Brown wrote:

> It looks like it is called for every quadrature point, probably the most
> expensive part of FE::reinit.

But look at which overloaded function is getting called:

> #2  0x00007f48c79e9fcf in FE<3u, (libMeshEnums::FEFamily)0>::reinit (this=0x252dbe0, elem=0x20f00a0, s=0, tolerance=9.9999999999999995e-07) at src/fe/fe_boundary.C:135

That's every quadrature point on reinit(elem, side), not on the
interior reinit(elem).  The inverse map gets used on sides to
transform a quadrature rule on a d-1 dimensional side element into
coordinates on the d-dimensional master element.  That may get
expensive for DG/jump calculations where every internal side is being
integrated on, but I'd be quite surprised if it was a major expense
for a CG assembly where the only side terms are on the domain
boundary.
---
Roy

On Thu, 11 Jun 2009, Jed Brown wrote:

>> That's every quadrature point on reinit(elem, side), not on the
>> interior reinit(elem).  The inverse map gets used on sides to
>> transform a quadrature rule on a d-1 dimensional side element into
>> coordinates on the d-dimensional master element.  That may get
>> expensive for DG/jump calculations where every internal side is being
>> integrated on, but I'd be quite surprised if it was a major expense
>> for a CG assembly where the only side terms are on the domain
>> boundary.
>
> Run ex4 and look at the profiling output.  An example run is attached.
> This shows 1 second total spent in element setup with almost half of
> that in inverse_map.  It's remarkably expensive, and called multiple
> times (not an even multiple) per boundary face.  That mesh has 6*32^2 =
> 6144 boundary faces, inverse_map is called more than 4 times this much
> (in one assembly).

It's called once per boundary quadrature point (note that the
vectorize inverse_map() just calls the one-point inverse_map() in a
loop) - that's 9 times per boundary face by default in ex4, I think,
which seems to square with your observations of ~4.5 local calls per
global boundary face when running on 2 processors.

> It still only amounts to 10% of assembly time, so this isn't an
> issue anywhere near the scale Tim is running into.

No.  Tim's also running some larger scale problems, where any
inefficiency on the boundary will scale up more slowly than the rest
of the assembly.

> Note that inverse_map could easily be a major bottleneck for CG in
> complex geometry, say bone structure analysis.

Possibly.  But it's been years since Ben decided that using
inverse_map for boundary quadrature was a good way to optimize for
programmer time at the expense of CPU time; IIRC Lorenzo just recently
became the first person to really disagree.

> In any case, I doubt this has anything to do with ghosted vectors.

No.
---
Roy

Dear Jed and Roy,

Please find attached now the result for 8 cores on 4 nodes.  Ghosted 
was still enabled.  I will now start the non-ghosted equivalent to 
this.

Roy: Did you have a look at the patch that I sent you yesterday?

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Fri, 12 Jun 2009, Tim Kroeger wrote:

> Roy: Did you have a look at the patch that I sent you yesterday?

Just got to it now.  I don't want to commit it until I can track down
what looks like a bug - on ex6 I got it to report that
build_inf_elem() takes 0.2663 seconds without subroutines included or
0.0000 seconds with them included.

But also, I'd like to wait and see if any other developers have strong
opinions.  This patch makes the PerfLog much more useful, but also
stretches the output to (typically) 107-110 characters.  Personally I
think we should just assume everyone has a 120+ char terminal, add
still another column for "Percent of Active Time With Sub" (hopefully
abbreviated better) and commit it.
---
Roy

On Fri, Jun 12, 2009 at 5:03 PM, Roy Stogner<roystgnr@...> wrote:
>
> On Fri, 12 Jun 2009, Tim Kroeger wrote:
>
>> Roy: Did you have a look at the patch that I sent you yesterday?
>
> Just got to it now.  I don't want to commit it until I can track down
> what looks like a bug - on ex6 I got it to report that
> build_inf_elem() takes 0.2663 seconds without subroutines included or
> 0.0000 seconds with them included.
>
> But also, I'd like to wait and see if any other developers have strong
> opinions.  This patch makes the PerfLog much more useful, but also
> stretches the output to (typically) 107-110 characters.  Personally I
> think we should just assume everyone has a 120+ char terminal, add
> still another column for "Percent of Active Time With Sub" (hopefully
> abbreviated better) and commit it.

That's fine with me.  There's no reason we need to stick to a certain
number of columns...

-- 
John

Agreed.

Derek

On Jun 12, 2009, at 4:10 PM, John Peterson wrote:

> On Fri, Jun 12, 2009 at 5:03 PM, Roy  
> Stogner<roystgnr@...> wrote:
>>
>> On Fri, 12 Jun 2009, Tim Kroeger wrote:
>>
>>> Roy: Did you have a look at the patch that I sent you yesterday?
>>
>> Just got to it now.  I don't want to commit it until I can track down
>> what looks like a bug - on ex6 I got it to report that
>> build_inf_elem() takes 0.2663 seconds without subroutines included or
>> 0.0000 seconds with them included.
>>
>> But also, I'd like to wait and see if any other developers have  
>> strong
>> opinions.  This patch makes the PerfLog much more useful, but also
>> stretches the output to (typically) 107-110 characters.  Personally I
>> think we should just assume everyone has a 120+ char terminal, add
>> still another column for "Percent of Active Time With Sub" (hopefully
>> abbreviated better) and commit it.
>
> That's fine with me.  There's no reason we need to stick to a certain
> number of columns...
>
> -- 
> John
>
> ------------------------------------------------------------------------------
> Crystal Reports - New Free Runtime and 30 Day Trial
> Check out the new simplified licensing option that enables unlimited
> royalty-free distribution of the report engine for externally facing
> server and web deployment.
> http://p.sf.net/sfu/businessobjects
> _______________________________________________
> Libmesh-devel mailing list
> Libmesh-devel@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-devel

Dear Jed and Roy,

Now the according result for non-ghosted vectors.

It seems as if nearly all of the performance difference is in the 
assembly methods.  Well, this is actually not very surprising, because 
the other big part of the application, that is solving the systems, 
does not deal with ghosted vectors at all.

The performance loss in the system assembly due to ghosted vectors is 
quite exactly the factor of 3 that I win back by the possibility to 
use 6 cores per node rather than 2.

Is there any idea what can be done?

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy,

On Fri, 12 Jun 2009, Roy Stogner wrote:

> On Fri, 12 Jun 2009, Tim Kroeger wrote:
>
>> Roy: Did you have a look at the patch that I sent you yesterday?
>
> Just got to it now.  I don't want to commit it until I can track down
> what looks like a bug - on ex6 I got it to report that
> build_inf_elem() takes 0.2663 seconds without subroutines included or
> 0.0000 seconds with them included.

Good point.

> Personally I
> think we should just assume everyone has a 120+ char terminal, add
> still another column for "Percent of Active Time With Sub" (hopefully
> abbreviated better) and commit it.

Good point.

I cared about both good points now, see attachment.

(The attachment also adds a log item in 
PetscVector::map_global_to_local_index(), that might help finding out 
whether my performance loss is really completely in that method.  You 
might or might not want to commit that as well.)

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:
> Dear Jed and Roy,
> 
> Now the according result for non-ghosted vectors.
> 
> It seems as if nearly all of the performance difference is in the
> assembly methods.  Well, this is actually not very surprising, because
> the other big part of the application, that is solving the systems, does
> not deal with ghosted vectors at all.

Yes, there is a bit more communication with ghosted vectors, which is
not normally needed, but that isn't using nearly enough time to
contribute.  It looks like all the time difference is in some unprofiled
events in assembly.  Are you sure there is nothing in your code that
branches based on whether ghosted vectors are being used?

In situations like these, I fall back to using conventional profiling
tools.  Unfortunately gprof can produce very misleading output, so I
almost never use it.  Since the time seems to not be in communication,
you could use callgrind (valgrind --tool=callgrind) and still have
relevant numbers when looking at the assembly function.  This option,
viewed with kcachegrind, is probably the most user-friendly thing and
should immediately show you where the time is spent.  You get source
annotation of the number of cycles spent on each line with a nice
graphical browser.  However, valgrind will be about 1000 times slower
than normal so you have to make the problem size truly tiny (and this
would disrupt things if there was e.g. an O(n^2) algorithm lurking
somewhere).

An alternative that requires no recompilation is Oprofile which you can
just start on one compute node, run the code, and look at the output.
Use opcontrol's --callgraph option.  This profiler is based on kernel
interrupts and has very little overhead (usually less than 5%).  Because
of the level it operates on, you need root to start the profiler.  You
can view the results in kcachegrind, but it's not as well integrated as
callgrind.  Regardless, you get very useful source annotation.

I think Eclipse has plugins for both valgrind and oprofile if you use
that.

Many clusters have commercial profiling tools as well, you may want to
ask the staff.  Helping track down these sort of performance bugs is
often part of their job description.  Sun Studio is free and has a
performance analyzer that is supposed to be pretty good, but I don't
have any direct experience with it.

Do yourself a favor and add a command-line option to stop the simulation
after N steps.  Do all your profiling with N<10 (a few minutes is plenty
of time to get accurate profiling).

Good luck.

Jed

Dear Jed,

On Mon, 15 Jun 2009, Jed Brown wrote:

> Tim Kroeger wrote:
>>
>> It seems as if nearly all of the performance difference is in the
>> assembly methods.  Well, this is actually not very surprising, because
>> the other big part of the application, that is solving the systems, does
>> not deal with ghosted vectors at all.
>
> In situations like these, I fall back to using conventional profiling
> tools.

Thank you for your detailed hints about profiling tools.  I'm somehow 
always very reluctant towards them for a number of reasons.  Skaling 
the whole application down such that it is faster by a factor of 1000 
could result in hard work -- and could also spoil the results.

Being root will be practically impossible.  The admin might be willing 
to help me, but he is located quite far away, so any help offered is 
based on email (or possible telephone), and he won't supply the root 
password to me.

Well, at least I have now enabled the already existing option of 
stopping the simulation after a few number of steps.  The problem is 
that even 10 steps already take about 5 hours.

I've added quite fine-grained START_LOG()/STOP_LOG() pairs now inside 
one of my assembly functions (that is, the one that is most probably 
responsible for the poor performance) and configured the application 
to run 4 steps only (which I suppose to take about two hours).  I'll 
keep you informed about the results.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Tim Kroeger wrote:

> Thank you for your detailed hints about profiling tools.  I'm somehow
> always very reluctant towards them for a number of reasons.  Skaling the
> whole application down such that it is faster by a factor of 1000 could
> result in hard work -- and could also spoil the results.

My experience is that this always pays off.  You aren't running the
small cases for the results it produces, you are running them to check
for code correctness.  My top three guidelines for developing parallel
code are

1. Have a single runtime parameter for the size of the problem, it
should scale between <1 second of runtime and a full production run.
(This can be an input mesh or some other parameter controlling problem
size.)

2. Make sure it works correctly in serial before trying in parallel.

3. Run in parallel on your workstation (small problem size) to check
correctness before moving to the cluster.


Also, I think this one is critical for any PDE solver:

* Manufacture solutions so that you have exact solutions to compare to.

This is really easy, even for very complex codes, if you write the
problem as F(u)=0 or F(u',u,t) = 0 and your code can accomodate
arbitrary forcing terms.  Choose the solution u(x,t) *before* choosing
the domain, boundary conditions or forcing.  The only requirement is
that it have sufficiently rich derivatives.  Products of transcendental
functions like tanh are good.  Then, using a symbolic algebra package
(Mathematica, Maple, Maxima, Sympy), apply your nonlinear differential
operators to manufacture a forcing term and print this as C code (these
packages can do this, don't worry if the expressions are pages long, you
never have to read them).  Paste the exact solutions and forcing terms
into your code and use the exact solutions for inhomogeneous boundary
conditions.  Now you can compare to highly nontrivial exact solutions.
Don't worry that they don't look anything like your real solutions
because the forcing terms are highly nonphysical, they *will* test that
your code is correct (converges to highly nontrivial exact solutions at
the correct rate).  This is way more useful than nearly degenerate
(physical) exact solutions that are commonly used for testing
correctness.  [/pulpit]

> Being root will be practically impossible.  The admin might be willing
> to help me, but he is located quite far away, so any help offered is
> based on email (or possible telephone), and he won't supply the root
> password to me.

It's definitely worth asking him what profiling tools are available.
Gprof is also an option.

> Well, at least I have now enabled the already existing option of
> stopping the simulation after a few number of steps.  The problem is
> that even 10 steps already take about 5 hours.

This doesn't make sense to me.  The full run was 576 steps in 18 hours
(ghosted case), which works out to about 2 minutes per step.  (Well,
that many assemblies, are there many assemblies per step?)

> I've added quite fine-grained START_LOG()/STOP_LOG() pairs now inside
> one of my assembly functions (that is, the one that is most probably
> responsible for the poor performance)

This is the poor-man's profiler.  It's a bit clumsy, but should work.

Jed

Tim Kroeger wrote:

> Well, there are in my case two or three parameters for this, but it's
> always again difficult to figure out where they have to be set, and it's
> always difficult not to forget to restore them once the profiling has
> been finished.

*runtime* parameters.  There is nothing to restore.

> I principally know this idea, but in my case it's really difficult,
> because I have fixed boundary values on a parts of a complicated
> geometry.  Parts of the geometry are hardcoded, and changing them for
> test purposes would mean to replace the code with a different one and
> could introduce additional bugs.  Also, external forcing terms are not
> supported in my application since they do not physically arise. Rather,
> things are driven by the boundary conditions (which are constant
> Dirichlet conditions on a complicated geometry), non-constant
> coefficients, and the fact that one of the equation gets a right hand
> side that is a function of the result of another equation.

It doesn't matter what the domain is.  You have to be able to
accommodate external forcing, but that is usually not so much work.  The
method is plenty general and can definitely be applied to your problem,
although I understand that it may be a nontrivial amount of work due to
legacy.  FWIW, I have no confidence in my own codes until this step has
been performed.  I'll stop trying to convince you now.

> In my assemble method, I have an element loop, and inside this a loop
> over the quadrature points.  main-loop-01-get-params is performed once
> for each quadrature point and will compute the value of the old solution
> at this point as well as the value of a couple of material parameters
> and other quantities (all of which depend on the solution of this or
> other systems but have been precomputed and stored in ExplicitSystem
> objects on their own).  This works mainly as follows:
> 
> Number val1 = 0.0;
> for (unsigned int l=0; l<phi_val1.size(); l++)
>   {
>     val1 +=
> phi_val1[l][qp]*(*system1.current_local_solution)(dof_indices_val1[l]);
>   }
> 
> There are about ten loops of this type.
[snip]
> Anyway, the main question is: Is it inefficient to call PETSc's
> VecGetArray()/VecRestoreArray() (and, in the ghosted case, also
> VecGhostGetLocalForm()/VecGhostRestoreLocalForm()) for each call to
> PetscVector::operator() if a large number of such calls are being done? 
> I guess it is not, and I suspect that doing a more efficient thing here
> could speed up my simulation considerably.  Jed, do you agree?

Yes, calling these in an inner loop is expensive.  The VecGhost
functions do RTTI (string comparison).  Calling them once per element
won't be a big deal, but it is definitely expensive when you use it for
every basis function at every quadrature point.  You should really just
get these arrays (as STL vectors if you like).

Jed

On Tue, 16 Jun 2009, Jed Brown wrote:

>> Anyway, the main question is: Is it inefficient to call PETSc's
>> VecGetArray()/VecRestoreArray() (and, in the ghosted case, also
>> VecGhostGetLocalForm()/VecGhostRestoreLocalForm()) for each call to
>> PetscVector::operator() if a large number of such calls are being done?
>> I guess it is not, and I suspect that doing a more efficient thing here
>> could speed up my simulation considerably.  Jed, do you agree?
>
> Yes, calling these in an inner loop is expensive.  The VecGhost
> functions do RTTI (string comparison).  Calling them once per element
> won't be a big deal, but it is definitely expensive when you use it for
> every basis function at every quadrature point.  You should really just
> get these arrays (as STL vectors if you like).

I think the PetscVector code is just being paranoid about guessing
what the get/restore "rules" are - knowing that we should call a
restore function "once one is finished using the object" or "when you
no longer need access to the array" seems unhelpful - *those*
conditions generally don't occur until the program is about to exit
(or at least until the vector is about to be destroyed) anyway.  It
wouldn't be hard to cache VecGetArray and VecGhostGetLocalForm values,
if we knew how long we can safely cache them for.  Does the cache need
to be renewed when VecAssembly/GhostUpdate Begin/End are called?  I
presume it wouldn't be invalidated by VecSetValues, VecAXPY, VecCopy,
etc.? 
---
Roy

Roy Stogner wrote:
> 
> On Tue, 16 Jun 2009, Jed Brown wrote:
> 
>>> Anyway, the main question is: Is it inefficient to call PETSc's
>>> VecGetArray()/VecRestoreArray() (and, in the ghosted case, also
>>> VecGhostGetLocalForm()/VecGhostRestoreLocalForm()) for each call to
>>> PetscVector::operator() if a large number of such calls are being done?
>>> I guess it is not, and I suspect that doing a more efficient thing here
>>> could speed up my simulation considerably.  Jed, do you agree?
>>
>> Yes, calling these in an inner loop is expensive.  The VecGhost
>> functions do RTTI (string comparison).  Calling them once per element
>> won't be a big deal, but it is definitely expensive when you use it for
>> every basis function at every quadrature point.  You should really just
>> get these arrays (as STL vectors if you like).
> 
> I think the PetscVector code is just being paranoid about guessing
> what the get/restore "rules" are - knowing that we should call a
> restore function "once one is finished using the object" or "when you
> no longer need access to the array" seems unhelpful - *those*
> conditions generally don't occur until the program is about to exit
> (or at least until the vector is about to be destroyed) anyway.  It
> wouldn't be hard to cache VecGetArray and VecGhostGetLocalForm values,
> if we knew how long we can safely cache them for.  Does the cache need
> to be renewed when VecAssembly/GhostUpdate Begin/End are called?  I
> presume it wouldn't be invalidated by VecSetValues, VecAXPY, VecCopy,
> etc.? ---

After modifying values, you absolutely MUST restore the arrays.  The
main reason for this is that some values (like norms) can be cached (it
makes the solvers much simpler to write without unnecessarily
recomputing norms).  With PETSc native vector types, the vectors are
backed by an array and the array doesn't move, but PETSc is extensible
and you cannot guarantee that any Vec you end up with is backed by an
array.  I could also imagine Vecs that are array-backed, but the array
is on a GC'd heap and could thus be relocated.  You should get the array
at the beginning of residual evaluation/matrix assembly/whatever and
restore it at the end.  Between getting and restoring the array, you
should not be calling Vec* functions (although it is technically not a
problem as long as you are not modifying values).  Following these
guidelines are normally not hard because you have to update ghost values
at these points anyway.

Jed

Dear Jed and Roy,

On Tue, 16 Jun 2009, Jed Brown wrote:

>> Anyway, the main question is: Is it inefficient to call PETSc's
>> VecGetArray()/VecRestoreArray() (and, in the ghosted case, also
>> VecGhostGetLocalForm()/VecGhostRestoreLocalForm()) for each call to
>> PetscVector::operator() if a large number of such calls are being done?
>> I guess it is not, and I suspect that doing a more efficient thing here
>> could speed up my simulation considerably.  Jed, do you agree?
>
> Yes, calling these in an inner loop is expensive.  The VecGhost
> functions do RTTI (string comparison).  Calling them once per element
> won't be a big deal, but it is definitely expensive when you use it for
> every basis function at every quadrature point.  You should really just
> get these arrays (as STL vectors if you like).

Hence, I have attached a patch for this purpose.  It includes a change 
in ex9 that demonstrates how to use of the new API.  I tested it, and 
it turned out to speed up ex9 by about 10% (one processor, devel mode, 
ghosted enabled).  Roy, what do you think?

I'll see how that speeds up my application, both with and without 
ghosted.  I expect it to speed up by much more than 10% since I have 
numerous loops of this type.  (This, of course, also increases the 
work to replace all the loops...)

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Thu, 18 Jun 2009, Tim Kroeger wrote:

> On Tue, 16 Jun 2009, Jed Brown wrote:
>
>> Yes, calling these in an inner loop is expensive.  The VecGhost
>> functions do RTTI (string comparison).  Calling them once per element
>> won't be a big deal, but it is definitely expensive when you use it for
>> every basis function at every quadrature point.  You should really just
>> get these arrays (as STL vectors if you like).
>
> Hence, I have attached a patch for this purpose.  It includes a change in ex9 
> that demonstrates how to use of the new API.  I tested it, and it turned out 
> to speed up ex9 by about 10% (one processor, devel mode, ghosted enabled). 
> Roy, what do you think?

I definitely like the results, but the API feels like "giving up".
Ideally I'd like to fix the speed problems with operator(); in theory
by doing that (and by adding a compile-time option to make basic
NumericVector calls inline instead of virtual; more on that later) we
ought to be able to get rid of the operator() overhead entirely
without requiring users to change their code.

If we can't do that, then your patch is certainly better than nothing,
but I'd like to wait before adding it - I don't want to encourage
users to change their code until we're certain it's necessary.
---
Roy

On Fri, 19 Jun 2009, Tim Kroeger wrote:

> I tried to do this now, see attached patch.

At first glance, the core of this patch looks exactly like what I had
in mind.  I'll want to spend some time testing (and ideally get some
reassurance from Jed) before making such a major change to such a
critical class, but this looks good.

> Also, I fixed something that looked like a bug in
> PetscVector::swap().

Ah - we weren't swapping the status flags!  That never bit me since I
always swapped vectors of the same type and in the same state, but it
could definitely have behaved horribly in future code.  Good catch!

I don't like making NumericVector::swap a non-pure-virtual, though,
since it's incomplete as is and has to be rederived in all subclasses.
I'd prefer to make failure-to-rederive into an obvious compile-time
error rather than a subtle run-time error, even if that means a couple
redundant lines of code here and there.  On the other hand, how often
do we add a whole new linear algebra package?  I guess it's fine as
is; we just need to remember to start calling it from the Trilinos and
Laspack subclasses too.

> Also, inside PetscVector::map_global_to_local_index(), I called
> VecGetOwnershipRange() direclty once (rather than calling it via an
> additional virtual function call twice), which makes a small, but
> observable speedup (I tested that separately).

Good idea.

> However, I currently don't see how you want to get rid of the virtual
> function call overhead.  I kept my new API for the while, as it is
> still somehow faster (although the differce is less pronounced).

Just an idea I ran by Ben recently over the phone; I hadn't mentioned
it on the list yet.

Although I don't want to get rid of the ability to do (limited) mixing
and matching of different linear algebra packages in the same code,
most users want to just pick one package and use it for every relevant
object.  We could move the current base classes to NumericVectorBase
and SparseMatrixBase, make NumericVector and SparseMatrix into
compile-time-selected typedefs of e.g. PetscVector and PetscMatrix,
and then declare all the leaf class functions non-virtual.  AFAIK this
still allows a C++ compiler to call those functions virtually from
pointers/references to the base class but also allows the compiler to
inline them when they're called from pointers/references to the leaf
class.

This idea will need testing to make sure that it works and it's worth
anything, naturally.
---
Roy

Tim Kroeger wrote:

> I assumed that functions like VecGetOwnershipRange(), VecGetSize(), and
> VecGetLocalSize() are allowed to be called between VecGetArray() and
> VecRestoreArray() (or, likewise, between VecGetLocalForm() and
> VecRestoreLocalForm()), although I'm not completely sure about that.
> Jed, can you comment on that?

These functions do not mutate any state and do not look at the values in
the array, so I can't think of any case where this would fail.  (I
haven't read the patch carefully, but the assumptions you mention above
are fine.)  Calling functions like VecScale or VecNorm while the
array/local form is gotten may not be safe.

Jed

On Fri, 19 Jun 2009, Roy Stogner wrote:

> On Fri, 19 Jun 2009, Tim Kroeger wrote:
>
>> I tried to do this now, see attached patch.
>
> At first glance, the core of this patch looks exactly like what I had
> in mind.  I'll want to spend some time testing (and ideally get some
> reassurance from Jed) before making such a major change to such a
> critical class, but this looks good.

Okay.  In any case, what will happen to the new method 
NumericVector::get() that is still in the patch?  Will you keep it or 
not?  I'm asking because I already use it in part of my application. 
If you will not include it in the library, I should revert the changes 
in my application as soon as possible (since in a couple of weeks I 
might have forgotten what I did exactly).  On the other hand, if you 
keep it, I see no reason for reverting my code since the new method is 
in any case faster (although, if you manage to avoid the virtual 
function call overhead, the difference will reduce to the check that 
_array_is_present is already true, which should really be negligible).

>> However, I currently don't see how you want to get rid of the virtual
>> function call overhead.  I kept my new API for the while, as it is
>> still somehow faster (although the differce is less pronounced).
>
> Just an idea I ran by Ben recently over the phone; I hadn't mentioned
> it on the list yet.
>
> Although I don't want to get rid of the ability to do (limited) mixing
> and matching of different linear algebra packages in the same code,
> most users want to just pick one package and use it for every relevant
> object.  We could move the current base classes to NumericVectorBase
> and SparseMatrixBase, make NumericVector and SparseMatrix into
> compile-time-selected typedefs of e.g. PetscVector and PetscMatrix,
> and then declare all the leaf class functions non-virtual.  AFAIK this
> still allows a C++ compiler to call those functions virtually from
> pointers/references to the base class but also allows the compiler to
> inline them when they're called from pointers/references to the leaf
> class.

I'm not quite sure about this.  In particular, I would think that 
inlining is only possible if you're calling from a local instance of 
the leaf class, but not from a pointer/reference to such a class. 
The reason is that I don't see any possibility for the compiler to see 
that your leaf class is really a leaf class, i.e. any user code could 
inherit further down and overload virtual functions and pass a 
pointer/refernce to such a class to the basic library.  C++ is missing 
a syntax for explicitly disallowing inheritance from a class, isn't 
it?

Even if it works theoretically, I see another problem: Of what type 
will, say, LinearImplicitSystem::solution be?  Will it be (an AutoPtr 
to) NumericVectorBase or NumericVector?  In either case, you can't 
have both the flexibility to mix solver algebra packaes and the 
possibility to inline virtual functions.

Anyway, I'd like to be proven wrong.  Also, compile-time selection of 
the linear algebra package would be fine for me since I'm only using 
PETSc anyway.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy/Jed,

I'm just running my complete application with the new patch installed, 
and it seems to me that it has become *slower*.  (It's not finished 
yet, but I can extrapolate the time.)

My question: Could it be that --enable-perflog and/or -log_summary 
reduce the performance essentially (that is, by 10% or more)?

Otherwise, things seem to me similar to those in George Orwell's 
"1984", where the weekly chocolate ration is *increased* from 30 grams 
to 25 grams.  (-:

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Mon, 22 Jun 2009, Tim Kroeger wrote:

> I'm just running my complete application with the new patch installed, and it 
> seems to me that it has become *slower*.  (It's not finished yet, but I can 
> extrapolate the time.)
>
> My question: Could it be that --enable-perflog and/or -log_summary reduce the 
> performance essentially (that is, by 10% or more)?

10% seems a little high, but yes, --enable-perflog reduces
performance; you're doing string matching every time you call it and
that can be expensive if you're calling it in a frequently used and
otherwise fast function.
---
Roy

On Mon, 22 Jun 2009, Tim Kroeger wrote:

> Okay.  In any case, what will happen to the new method NumericVector::get() 
> that is still in the patch?  Will you keep it or not?

I'll probably leave out the ex3 usage, or maybe move that to a more
advanced example, but there's no reason not to keep the method in the
library itself.  It's noticeably more efficient in every use case
right now and it still would be in some use cases in the future even if
my wacky idea worked.

> I'm not quite sure about this.  In particular, I would think that inlining is 
> only possible if you're calling from a local instance of the leaf class, but 
> not from a pointer/reference to such a class. The reason is that I don't see 
> any possibility for the compiler to see that your leaf class is really a leaf 
> class, i.e. any user code could inherit further down and overload virtual 
> functions and pass a pointer/refernce to such a class to the basic library. 
> C++ is missing a syntax for explicitly disallowing inheritance from a class, 
> isn't it?

Yeah, it looks like you're right.  I'd assumed that removing the
"virtual" keyword in a leaf class would be enough to tell the compiler
to use non-virtual dispatch when calling that method through a pointer
to the leaf class, but when actually testing that out, it doesn't work
that way.

Of course, we could get the same behavior by making a configure-time
option that in addition to renaming classes turns some of their
virtual keywords on or off as appropriate, but that's much more work
than what I'd originally planned, so it'll go on the back burner.

> Even if it works theoretically, I see another problem: Of what type will, 
> say, LinearImplicitSystem::solution be?  Will it be (an AutoPtr to) 
> NumericVectorBase or NumericVector?  In either case, you can't have both the 
> flexibility to mix solver algebra packaes and the possibility to inline 
> virtual functions.

Not at the same time; you'd have to choose one or the other from
./configure.
---
Roy

Dear Roy and Jed,

Now, I have performed a couple of test runs of my applications with 
the new implementation of PetscVector, both with and without ghosted 
vectors.  All runs were performed in devel mode and without 
performance logging (to make it comparable to the runs that were 
performed earlier).  I prepared a list for you here of these runs as 
well as the earliear runs (with old PetscVector implementation):

PetVec	ghost	#cores	#nodes	hh:mm:ss
----------------------------------------
old	no	8	3	11:34:10
old	no	8	3	11:35:54
old	ghost	8	3	17:25:28
old	ghost	8	3	16:33:23
old	ghost	24	4	10:42:34
old	ghost	24	4	11:05:30
old	ghost	24	4	11:24:17
new	no	8	3	12:10:44
new	no	9	3	12:01:35
new	ghost	8	3	11:59:33
new	ghost	9	3	11:21:11
new	ghost	18	3	08:31:03
new	ghost	24	4	08:37:35

I also looked at the results, and they coincide up to the usual 
numerical differences.

It seems as if the new PetscVector implementation makes the 
non-ghosted code slower.  I have no explanation for this.  On the 
other hand, the ghosted code is much faster now than before and 
compares well with the non-ghosted plus old PetscVector.  And for 
applications where the number of employed cores per node is limited 
due to memory requirements and where a major part of the time is spent 
in system assembly (which both is true for my application), an 
essential speedup is possible now using ghosted plus new PetscVector 
(and thus more cores per node).

I would vote for making both default soon.

Any opinions?

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Mon, 29 Jun 2009, Tim Kroeger wrote:

> Now, I have performed a couple of test runs of my applications with the new 
> implementation of PetscVector, both with and without ghosted vectors.  All 
> runs were performed in devel mode and without performance logging (to make it 
> comparable to the runs that were performed earlier).  I prepared a list for 
> you here of these runs as well as the earliear runs (with old PetscVector 
> implementation):

Sorry to get back to you so late on this.  I think my subconscious was
stalling, because I can still hardly write the necessary reply without
wincing:

Can you run all those again, in opt mode?

METHOD=devel is intended to be a compromise between the "libmesh
assertions are disabled and debugging the code takes forever" behavior
of METHOD=opt and the "g++ STL checking is enabled and running the
code takes forever" behavior of METHOD=dbg.  It's great for catching
issues that trip a libmesh_assert() before they segfault in an out of
bounds vector access.  But it's probably worthless for performance
timing.  It's not nearly as slow as dbg mode, but it's significantly
slower than normal and the degree of slowdown depends in part on just
how paranoid the asserts in a particular code path are.

> It seems as if the new PetscVector implementation makes the non-ghosted code 
> slower.  I have no explanation for this.

I'll bet this behavior goes away with METHOD=opt.  If so, then I'd
agree that we should make the ghosted code the default.
---
Roy

Dear Roy,

On Thu, 2 Jul 2009, Roy Stogner wrote:

> Sorry to get back to you so late on this.  I think my subconscious was
> stalling, because I can still hardly write the necessary reply without
> wincing:
>
> Can you run all those again, in opt mode?

Seems as if you are right.  I will do this, but it will of course take 
some time.  Luckily enough, I didn't carry on working on my 
application in the meantime.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

On Fri, 3 Jul 2009, Tim Kroeger wrote:

> Dear Roy,
>
> On Thu, 2 Jul 2009, Roy Stogner wrote:
>
>> Sorry to get back to you so late on this.  I think my subconscious was
>> stalling, because I can still hardly write the necessary reply without
>> wincing:
>> 
>> Can you run all those again, in opt mode?
>
> Seems as if you are right.  I will do this, but it will of course take some 
> time.  Luckily enough, I didn't carry on working on my application in the 
> meantime.

Thanks.

Oh, and just to make sure: you'll need to build PETSc in optimized
mode too to get reasonable timing results.  Just this afternoon I
discovered that my code has been running very slowly on one system
because I'd been mistakenly defaulting to a debug mode PETSc module.
---
Roy

Dear Roy,

On Fri, 3 Jul 2009, Roy Stogner wrote:

> On Fri, 3 Jul 2009, Tim Kroeger wrote:
>
>> Seems as if you are right.  I will do this, but it will of course take some 
>> time.  Luckily enough, I didn't carry on working on my application in the 
>> meantime.

First result: Ghosted enabled, new PetscVector implementation, 18 
cores on 3 nodes, opt mode: 8:03:03.  (Was 8:31:03 in devel mode.)

> Oh, and just to make sure: you'll need to build PETSc in optimized
> mode too to get reasonable timing results.  Just this afternoon I
> discovered that my code has been running very slowly on one system
> because I'd been mistakenly defaulting to a debug mode PETSc module.

Good that you mention it.  Of course I had forgotten that.  Actually, 
I have never so far used PETSc in optimized mode.  My application is 
far from being in a final state, and during development the debug 
version of PETSc is somehow useful for backtracking crashes.  Anyway, 
for sensible timing results, you're probably right that PETSc's 
optimized mode should be prefered (although, as we have seen, my 
application does actually not spend too much time inside PETSc at 
all).

Anything else that I might have forgotten?

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany

Dear Roy,

Now, the more-or-less-ultimate results.  Both PETSc and libMesh were 
compiled in optimized mode.

PetVec	ghost	#cores	#nodes	hh:mm:ss
----------------------------------------
new	ghost	18	3	07:02:54
new	ghost	9	3	08:19:55
new	no	9	3	08:35:28
old	no	9	3	08:21:52
old	ghost	9	3	crashed
old	ghost	18	3	crashed

Hence, for ghosted disabled, new and old PetscVector implementations 
are about equally fast.  Also, with the new PetscVector 
implementation, enabling ghosted on the same number of cores and nodes 
keeps the performance comparable, and increasing the number of cores 
speeds things up essentially.  Strange is the obervation that the old 
implementation with ghosted enabled now crashes my application.  The 
crash occurs in the post-processing step at the end (which I had 
temporarily disabled, but I now re-enabled it).

I'll try to track down what happened.

Best Regards,

Tim

-- 
Dr. Tim Kroeger
tim.kroeger@...            Phone +49-421-218-7710
tim.kroeger@...            Fax   +49-421-218-4236

Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany