libMesh: A C++ Finite Element Library

On Tue, Oct 7, 2008 at 10:32 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear libMesh team,
>
> As far as I see, there seems not to be any interface in libMesh to
> solve a system in a matrix free way.
>
> I have to solve a linear system whose matrix basically has the
> following structure:
>
> A = B + v*w^T
>
> where B is a sparse matrix and v and w are vectors, so that v*w^T is a
> full matrix.  Obviously, it is not advisable to store the full matrix.
> Instead, I would like to store B, v, and w, and then supply a method
> to multiply the matrix with a given vector (which is obviously quite
> easy, even in parallel).
>
> What whould be the things to do?  I guess, first I have to find out
> PETSc's API for matrix free solving (provided that such a thing
> exists, but I think it does).  Then, I can either use PETSc directly
> or -- presumably better -- extend the API of your LinearSolver class
> towards matrix free solving.  Of course, the latter would require to
> do something for the other solver packages as well (Laspack etc.), but
> I would just call libmesh_not_implemented() there.
>
> Would such an API extension be welcome?  Do you have any different
> ideas or suggestions?

Definitely.  I'm not sure what your intended use of this functionality
will be, but solving "bordered systems" efficiently is one possible
application.  I'm not sure if we should call this "matrix free" since
I think that has other connotations for our developers in the context
of Jacobian free solutions of nonlinear systems.

One way to implement the desired functionality in PETSc is to define a
custom matrix-vector product operation.  See, for example,
MatShellSetOperation.  If at first you can only get this working for
PETSc, that's fine.  I'm not sure what the best way of abstracting
this concept is.

-- 
John

On Tue, 7 Oct 2008, John Peterson wrote:

> On Tue, Oct 7, 2008 at 10:32 AM, Tim Kroeger
> <tim.kroeger@...> wrote:
>> As far as I see, there seems not to be any interface in libMesh to
>> solve a system in a matrix free way.
>>
>> I have to solve a linear system whose matrix basically has the
>> following structure:
>>
>> A = B + v*w^T
>>
>> where B is a sparse matrix and v and w are vectors, so that v*w^T is a
>> full matrix.  Obviously, it is not advisable to store the full matrix.
>> Instead, I would like to store B, v, and w, and then supply a method
>> to multiply the matrix with a given vector (which is obviously quite
>> easy, even in parallel).
>>
>> What whould be the things to do?  I guess, first I have to find out
>> PETSc's API for matrix free solving (provided that such a thing
>> exists, but I think it does).

You can basically do a MatCreateShell and then hand it to any of their
existing APIs, right?  Just because PETSc isn't C++ doesn't mean it
isn't object oriented.  But I've never tried it myself, and I'll bet
your preconditioning options get tricky.

>> Then, I can either use PETSc directly
>> or -- presumably better -- extend the API of your LinearSolver class
>> towards matrix free solving.  Of course, the latter would require to
>> do something for the other solver packages as well (Laspack etc.), but
>> I would just call libmesh_not_implemented() there.
>>
>> Would such an API extension be welcome?  Do you have any different
>> ideas or suggestions?
>
> Definitely.

I agree.

> I'm not sure if we should call this "matrix free" since I think that
> has other connotations for our developers in the context of Jacobian
> free solutions of nonlinear systems.

This is also true.  I'm not sure if "shell matrix" is a standard term
or a PETSc-ism, but I'd be happy to adopt it.

> One way to implement the desired functionality in PETSc is to define a
> custom matrix-vector product operation.  See, for example,
> MatShellSetOperation.  If at first you can only get this working for
> PETSc, that's fine.  I'm not sure what the best way of abstracting
> this concept is.

Me neither.  My first gut reaction is to create a ShellMatrix subclass
of SparseMatrix.
---
Roy

There is a way to do matrix free solving with libMesh and Petsc.  Take  
a look at the PetscNonlinearSolver (which you can get by using a  
NonlinearImplicitSystem).  If you just set the residual() function and  
not the jacobian() then you are solving matrix free using Petsc SNES.

Note that there is a "matvec" function on PetscNonlinearSolver  
(inherited from NonlinearSolver).  I'm not quite sure how to use it  
though... it's something Ben put in I believe.

If you want to turn on that "matvec" function in PetscNonlinearSolver  
and make it into a MatShell or PCShell that might be a good way to go.

Derek

On Oct 7, 2008, at 9:32 AM, Tim Kroeger wrote:

> Dear libMesh team,
>
> As far as I see, there seems not to be any interface in libMesh to
> solve a system in a matrix free way.
>
> I have to solve a linear system whose matrix basically has the
> following structure:
>
> A = B + v*w^T
>
> where B is a sparse matrix and v and w are vectors, so that v*w^T is a
> full matrix.  Obviously, it is not advisable to store the full matrix.
> Instead, I would like to store B, v, and w, and then supply a method
> to multiply the matrix with a given vector (which is obviously quite
> easy, even in parallel).
>
> What whould be the things to do?  I guess, first I have to find out
> PETSc's API for matrix free solving (provided that such a thing
> exists, but I think it does).  Then, I can either use PETSc directly
> or -- presumably better -- extend the API of your LinearSolver class
> towards matrix free solving.  Of course, the latter would require to
> do something for the other solver packages as well (Laspack etc.), but
> I would just call libmesh_not_implemented() there.
>
> Would such an API extension be welcome?  Do you have any different
> ideas or suggestions?
>
> Best Regards,
>
> Tim
>
> -- 
> Dr. Tim Kroeger                                        Phone  
> +49-421-218-7710
> tim.kroeger@..., tim.kroeger@...  Fax    
> +49-421-218-4236
>
> MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany
>
> Amtsgericht Bremen HRB 16222
> Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen
>
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's  
> challenge
> Build the coolest Linux based applications with Moblin SDK & win  
> great prizes
> Grand prize is a trip for two to an Open Source event anywhere in  
> the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users@...
> https://lists.sourceforge.net/lists/listinfo/libmesh-users

Dear all,

thank you very much for all your suggestions.

I started implementing the thing now.  Find attached a patch that 
contains my first attempt.  It compiles well at the moment, but it 
most probably will not work as expected, because I did not have the 
time to test it yet, and I am not too much experienced in PETSc 
programming.

However, I would appreciate your feedback concerning my API 
suggestion.  Let me explain this to you:

Coarse scale: LinearImplicitSystem gets a method attach_shell_matrix() 
that the user should call if he wants to do "matrix free" linear 
solving.  What then happens is that LinearImplicitSystem::solve() 
notices that there is a shell matrix and then calls a different solver 
function.

Fine scale: LinearSolver gets another solve() method that eats a shell 
matrix function rather than a SparseMatrix object.  This calls the 
respective PETSc functions to do a matrix free solve.

It is a little bit tricky not to get confused with the contexts that 
these shell functions need: The user supplied shell function has a 
reference to EquationSystems; that should be enough as a context.  The 
shell function that LinearImplicitSystem passes to LinearSolver cannot 
do this because (a) LinearSolver does not know about EquationSystems, 
and (b) it wouldn't know which system is the right one.  Hence, I use 
a pointer to *this here, casted to void*.  The shell function that is 
passed to PETSc needs both the user shell function and user context 
(which is the pointer to the LinearImplicitSystem of course), hence I 
construct a std::vector<void*>* which I again cast into void*.  I 
agree that this is somehow confusing, but I think it should work.

Any comments are welcome.  I will test the thing in the next days (by 
implementing the application that I need it for) and let you know when 
it works.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Wed, Oct 8, 2008 at 10:52 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear all,
>
> thank you very much for all your suggestions.
>
> I started implementing the thing now.  Find attached a patch that contains
> my first attempt.  It compiles well at the moment, but it most probably will
> not work as expected, because I did not have the time to test it yet, and I
> am not too much experienced in PETSc programming.

OK, I'd be interested to see the results of your testing.

> Any comments are welcome.  I will test the thing in the next days (by
> implementing the application that I need it for) and let you know when it
> works.

I'm a little nervous about casting a function pointer to a void*.
Technically this is not allowed (something about void* being reserved
for data pointers, google for "function pointer to void*" on
comp.lang.c++.moderated) which I assume is why you have used the
implementation-defined reinterpret_cast.  I think the main problem is
that "a void* is not required to be of adequate size to hold a
function pointer" but it is on all POSIX systems.

One possiblity might be instead of a vector<void*> for inner_ctx, a
pair<void*, void(*)(void)> unless you anticipate adding more and more
stuff to the inner_ctx vector...

This may be one of those things that is not allowed, but is
nevertheless done by everyone anyway?  We may have to do slightly
non-portable stuff to get around PETSc function prototype
requirements.  At the very least it might be worth checking the
results of the void* casts against NULL before you try to use them,
though I think only dynamic_cast is required to return NULL upon
failure.  Also, we might ask you to make this code
backwards-compatible with older versions of PETSc (or at least throw
errors) before we check it in.

-- 
John

On Wed, 8 Oct 2008, John Peterson wrote:

> I'm a little nervous about casting a function pointer to a void*.
> Technically this is not allowed (something about void* being reserved
> for data pointers, google for "function pointer to void*" on
> comp.lang.c++.moderated) which I assume is why you have used the
> implementation-defined reinterpret_cast.  I think the main problem is
> that "a void* is not required to be of adequate size to hold a
> function pointer" but it is on all POSIX systems.

Good catch.

> One possiblity might be instead of a vector<void*> for inner_ctx, a
> pair<void*, void(*)(void)> unless you anticipate adding more and more
> stuff to the inner_ctx vector...

How about a functor abstract base class?  Anyone?  We are using C++
still, right?  ;-)
---
Roy

On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner <roystgnr@...> wrote:
>
> On Wed, 8 Oct 2008, John Peterson wrote:
>
>> I'm a little nervous about casting a function pointer to a void*.
>> Technically this is not allowed (something about void* being reserved
>> for data pointers, google for "function pointer to void*" on
>> comp.lang.c++.moderated) which I assume is why you have used the
>> implementation-defined reinterpret_cast.  I think the main problem is
>> that "a void* is not required to be of adequate size to hold a
>> function pointer" but it is on all POSIX systems.
>
> Good catch.
>
>> One possiblity might be instead of a vector<void*> for inner_ctx, a
>> pair<void*, void(*)(void)> unless you anticipate adding more and more
>> stuff to the inner_ctx vector...
>
> How about a functor abstract base class?  Anyone?  We are using C++
> still, right?  ;-)

I briefly considered this but MatShellSetOperation does required a
pointer to a void(void) function to be passed.

http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatShellSetOperation.html

With a function object or even just a non-static member function of
the class, can this be done with some mem_fun tricks, e.g.

mem_fun(&Foo::the_function)

?

-- 
John

On Wed, 8 Oct 2008, John Peterson wrote:

> On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner <roystgnr@...> wrote:
>>
>> How about a functor abstract base class?  Anyone?  We are using C++
>> still, right?  ;-)
>
> I briefly considered this but MatShellSetOperation does required a
> pointer to a void(void) function to be passed.
>
> http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Mat/MatShellSetOperation.html

Oh, yeah, I wouldn't expect PETSc to take a C++ class.  We'd have some
behind-the-scenes trickery to make sure that PETSc got a C function
that ended up calling the correct C++ functor.

> With a function object or even just a non-static member function of
> the class, can this be done with some mem_fun tricks, e.g.
>
> mem_fun(&Foo::the_function)

If the C function takes a void*, everything's great - you have a
global function which casts that argument to a pointer to your
functor, and you call your functor.

If the C function really took void, things would be nasty (since you'd
need a separate C function for every shell matrix functor, but the
former have to be defined at compile time and the latter ought to be
createable at runtime), but it looks like the operator function ends
up getting cast to something which takes the same arguments (e.g.
Mat,Vec,Vec) as the corresponding PETSc matrix operation, and so we
could use MatShellGetContext to grab a pointer to our functor.
---
Roy

On Wed, Oct 8, 2008 at 12:55 PM, Jed Brown <jed@...> wrote:
> On Wed 2008-10-08 12:42, John Peterson wrote:
>> Thanks Jed...
>>
>> On Wed, Oct 8, 2008 at 12:24 PM, Jed Brown <jed@...> wrote:
>> >
>> > class MyShell {
>> > public:
>> >      int mult(...);
>> > };
>> >
>> > PetscErrorCode MyMult(Mat A,Vec x,Vec y)
>> > {
>> >  MyShell *shell;
>> >  MatShellGetContext(A,(void**)&shell);
>> >  shell->mult(x,y);
>> > }
>> >
>> >
>> > // elsewhere, perhaps in the MyShell constructor
>> > MatShellCreate(comm,m,n,M,N,&B);
>> > MatShellSetOperation(B,MATOP_MULT,(void(*)(void))MyMult);
>>
>>
>> MatShellCreate (or CreateShell?) takes a void* ctx, so if it were in
>> MyShell constructor we'd also pass static_cast<void*>(this) ?
>
> Yeah, that's what I meant.
>
> The same procedure is recommended for other callbacks (convergence
> tests, custom preconditioners, function/Jacobian evaluation, etc) in
> PETSc.  Presumably you do something like this with the SNES interface?

Actually, it looks like (src/numerics/petsc_nonlinear_solver.C) we use
some specially-named global functions for e.g. SNESSetJacobian, but we
do pass "this" as the context pointer, where "this" is LibMesh's
PetscNonlinearSolver class.    I'd much rather encapsulate the global
functions in a class as you have done.

-- 
John

On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner <roystgnr@...> wrote:
>
> On Wed, 8 Oct 2008, John Peterson wrote:
>
>> I'm a little nervous about casting a function pointer to a void*.
>> Technically this is not allowed (something about void* being reserved
>> for data pointers, google for "function pointer to void*" on
>> comp.lang.c++.moderated) which I assume is why you have used the
>> implementation-defined reinterpret_cast.  I think the main problem is
>> that "a void* is not required to be of adequate size to hold a
>> function pointer" but it is on all POSIX systems.
>
> Good catch.

At the risk of being overly pedantic, I've just added a configure test
which checks for the size of a function pointer.  I'm not convinced
that you will still need to do a cast of this nature, Tim, but if you
do, could you put throw something like the following:

#if LIBMESH_SIZEOF_FUNCTION_POINTER != LIBMESH_SIZEOF_VOID_P
std::cerr << "Appropriate error message." << std::endl;
libmesh_error();
#endif

in the code somewhere nearby?  It just feels like a nasty bug waiting
to happen otherwise...

-- 
John

Dear Roy and John and all,

Please find attached my next version.

I enclosed all my PETSc calles with "#if 
PETSC_VERSION_LESS_THAN(2,3,3)" and added error messages in the 
"#else" case since I am not familiar with PETSc's API history. 
Whenever somebody might need the shell matrix functionality together 
with older PETSc versions, he might want to implement that himself.

I switched the user API for shell matrix usage from a callback 
function to a class, which in fact seems to be much more sensible.  I 
made an abstract base class "ShellMatrix" and three derived classes. 
One of the derived classes just wraps a SparseMatrix.  I noticed that 
SparseMatrix apparently did not have a method to compute a 
matrix-vector product, so I added this functionality to SparseMatrix 
and PetscMatrix (whereas calling libmesh_not_implemented() in 
LaspackMatrix).

Thank you for the hint that a function pointer may not be casted into 
a void*.  I did not know this before.  Anyway, this is no longer 
necessary since I am now working with a pointer to a class.

Again, everything compiles well now but is not yet tested.  Comments 
are welcome.  I'll keep you informed about my progress.

By the way, why is the .depend file contained in the repository?

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Thu, Oct 9, 2008 at 7:04 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear Roy and John and all,
>
> Please find attached my next version.
>
> I enclosed all my PETSc calles with "#if PETSC_VERSION_LESS_THAN(2,3,3)" and
> added error messages in the "#else" case since I am not familiar with
> PETSc's API history. Whenever somebody might need the shell matrix
> functionality together with older PETSc versions, he might want to implement
> that himself.

OK, that may unnecessarily break your code when Petsc 2.3.4 (or
whatever their next version will be) comes out though.  I'd do
instead:

#if PETSC_VERSION_LESS_THAN(2,3,2)
// suitable error or unimplemented message
#else
// your current code
#endif


> I switched the user API for shell matrix usage from a callback function to a
> class, which in fact seems to be much more sensible.  I made an abstract
> base class "ShellMatrix" and three derived classes. One of the derived
> classes just wraps a SparseMatrix.  I noticed that SparseMatrix apparently
> did not have a method to compute a matrix-vector product, so I added this
> functionality to SparseMatrix and PetscMatrix (whereas calling
> libmesh_not_implemented() in LaspackMatrix).

This functionality (matrix-vector product) is actually in the
NumericVector class, it's the add_vector routine which takes a
SparseMatrix.  Enough people have been confused about this that we
should probably implement it in SparseMatrix as well, but I'd vote for
an implementation in terms of the existing NumericVector one (if
possible) rather than duplicating code.


> By the way, why is the .depend file contained in the repository?

A couple of reasons: this way we don't require all users to have perl
installed (though I don't know how they wouldn't...) to generate
dependencies.  Second, we used to generate dependencies with the
compiler instead of perl and it really took forever, so it made sense
to check in the .depend file.

-- 
John

John Peterson wrote:
> On Thu, Oct 9, 2008 at 7:04 AM, Tim Kroeger
> <tim.kroeger@...> wrote:
>> By the way, why is the .depend file contained in the repository?
> 
> A couple of reasons: this way we don't require all users to have perl
> installed (though I don't know how they wouldn't...) to generate
> dependencies.  Second, we used to generate dependencies with the
> compiler instead of perl and it really took forever, so it made sense
> to check in the .depend file.

A third reason that makes sense to me: when someone does an "svn 
update", the dependencies can change.  Either that update needs to make 
the .depend file change with them, or the .depend file can't be trusted 
and we need to regenerate it on every build.

I'd say we're open to better ideas, though.  In particular, I've just 
noticed that our example (and my application) .depend files end up with 
the libMesh path hardcoded into them, and so when switching an 
application from one LIBMESH_DIR to another you need to wipe and 
regenerate (the application's) .depend *anyway*.  If that gotcha ever 
catches me, I'll be tempted to either change around make_dependencies.pl 
or find a way to replace it entirely.
---
Roy

Dear John,

On Thu, 9 Oct 2008, John Peterson wrote:

>> I enclosed all my PETSc calles with "#if PETSC_VERSION_LESS_THAN(2,3,3)" and
>> added error messages in the "#else" case since I am not familiar with
>> PETSc's API history. Whenever somebody might need the shell matrix
>> functionality together with older PETSc versions, he might want to implement
>> that himself.
>
> OK, that may unnecessarily break your code when Petsc 2.3.4 (or
> whatever their next version will be) comes out though.  I'd do
> instead:
>
> #if PETSC_VERSION_LESS_THAN(2,3,2)
> // suitable error or unimplemented message
> #else
> // your current code
> #endif

Well, that's actually almost what I did, except that I compared to 
(2,3,3) rather than (2,3,2), since I concluded from other usages of 
that macro that it is a "<" comparison rather than a "<=".  (I agree 
that my mail could lead to the assumption that I had done it the wrong 
way around.)

>> I switched the user API for shell matrix usage from a callback function to a
>> class, which in fact seems to be much more sensible.  I made an abstract
>> base class "ShellMatrix" and three derived classes. One of the derived
>> classes just wraps a SparseMatrix.  I noticed that SparseMatrix apparently
>> did not have a method to compute a matrix-vector product, so I added this
>> functionality to SparseMatrix and PetscMatrix (whereas calling
>> libmesh_not_implemented() in LaspackMatrix).
>
> This functionality (matrix-vector product) is actually in the
> NumericVector class, it's the add_vector routine which takes a
> SparseMatrix.  Enough people have been confused about this that we
> should probably implement it in SparseMatrix as well, but I'd vote for
> an implementation in terms of the existing NumericVector one (if
> possible) rather than duplicating code.

Okay, done.  Also, I did the other way round for ShellMatrix, i.e. I 
added NumericVector::add_vector(ShellMatrix,NumericVector) which is 
implemented in terms of ShellMatrix::vector_mult_add().

My application that tests the shell matrix is still under 
construction.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

Dear John/Roy/all,

As Roy predicted, preconditioning is tricky.  I extended the 
ShellMatrix class (and all the derived classes) by a get_diagonal() 
method, thus allowing at least JACOBI preconditioning, which for the 
moment seems to be sufficient for my example.  This required to add an 
analogous function to the SparseMatrix class and a pointwise_mult() 
function to NumericVector.

A patch with the current version is attached.

It works for my example on one processor, where "works" means "does 
not crash and gives results that might or might not be correct".

I am currently unable to test it in parallel since there is some 
trouble with the cluster I've been using.  I hope this will be 
resolved soon.

I'll keep you informed.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Thu, 16 Oct 2008, Tim Kroeger wrote:

> As Roy predicted, preconditioning is tricky.  I extended the ShellMatrix 
> class (and all the derived classes) by a get_diagonal() method, thus allowing 
> at least JACOBI preconditioning, which for the moment seems to be sufficient 
> for my example.  This required to add an analogous function to the 
> SparseMatrix class and a pointwise_mult() function to NumericVector.

I'm a bit unhappy with the current inheritance tree - there's no
inheritance between ShellMatrix and SparseMatrix, even though the
latter is already practically an abstract base class that declares the
functions the former will need to use?  There's no getting around the
need to duplicate a lot of function definitions (get_diagonal() is
going to be different in the different classes no matter how things
are factored) but it seems like a little refactoring (make both
ShellMatrix and SparseMatrix inherit from a MatrixBase?) would at
least give us more concise code using them.

Feel free to ignore my vague opinions, though; the patch looks like a
great start and the shell matrix feature is a good idea regardless.
---
Roy

On Thu, Oct 16, 2008 at 10:20 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear John/Roy/all,
>
> As Roy predicted, preconditioning is tricky.  I extended the ShellMatrix
> class (and all the derived classes) by a get_diagonal() method, thus
> allowing at least JACOBI preconditioning, which for the moment seems to be
> sufficient for my example.  This required to add an analogous function to
> the SparseMatrix class and a pointwise_mult() function to NumericVector.
>
> A patch with the current version is attached.

Looks like you're learning more about PETSc than you ever wanted to
know :-)  The pointwise multiplication and diagonal getter
functionality is nice to have, thanks for contributing it back to the
library.

It's a little silly we have to do things like:

const_cast<void*>(static_cast<const void*>(&shell_matrix)),

in order to make use of the (supposedly better) C++ casts when the C-style

(void*)(&shell_matrix)

works and looks so much neater.  My preference is to stick with C++
casts since it does make it very explicit, but sheesh...


> It works for my example on one processor, where "works" means "does not
> crash and gives results that might or might not be correct".
>
> I am currently unable to test it in parallel since there is some trouble
> with the cluster I've been using.  I hope this will be resolved soon.
>
> I'll keep you informed.

OK, sounds good.

-- 
John

On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner <roystgnr@...> wrote:
>
> On Thu, 16 Oct 2008, Tim Kroeger wrote:
>
>> As Roy predicted, preconditioning is tricky.  I extended the ShellMatrix
>> class (and all the derived classes) by a get_diagonal() method, thus allowing
>> at least JACOBI preconditioning, which for the moment seems to be sufficient
>> for my example.  This required to add an analogous function to the
>> SparseMatrix class and a pointwise_mult() function to NumericVector.
>
> I'm a bit unhappy with the current inheritance tree - there's no
> inheritance between ShellMatrix and SparseMatrix, even though the
> latter is already practically an abstract base class that declares the
> functions the former will need to use?  There's no getting around the
> need to duplicate a lot of function definitions (get_diagonal() is
> going to be different in the different classes no matter how things
> are factored) but it seems like a little refactoring (make both
> ShellMatrix and SparseMatrix inherit from a MatrixBase?) would at
> least give us more concise code using them.

I think I'll withhold judgment until I've seen the new ShellMatrix
header file... It doesn't seem like they should necessarily be treated
as members of the SparseMatrix inheritance family.  It's bad enough
that we can't simultaneously keep the LaspackMatrix interface
up-to-date with changes to the Petsc one, and now we want to add one
more family member? ;-)

-- 
John

On Thu, Oct 16, 2008 at 11:20 AM, John Peterson <jwpeterson@...> wrote:
> On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner <roystgnr@...> wrote:
>>
>> On Thu, 16 Oct 2008, Tim Kroeger wrote:
>>
>>> As Roy predicted, preconditioning is tricky.  I extended the ShellMatrix
>>> class (and all the derived classes) by a get_diagonal() method, thus allowing
>>> at least JACOBI preconditioning, which for the moment seems to be sufficient
>>> for my example.  This required to add an analogous function to the
>>> SparseMatrix class and a pointwise_mult() function to NumericVector.
>>
>> I'm a bit unhappy with the current inheritance tree - there's no
>> inheritance between ShellMatrix and SparseMatrix, even though the
>> latter is already practically an abstract base class that declares the
>> functions the former will need to use?  There's no getting around the
>> need to duplicate a lot of function definitions (get_diagonal() is
>> going to be different in the different classes no matter how things
>> are factored) but it seems like a little refactoring (make both
>> ShellMatrix and SparseMatrix inherit from a MatrixBase?) would at
>> least give us more concise code using them.
>
> I think I'll withhold judgment until I've seen the new ShellMatrix
> header file... It doesn't seem like they should necessarily be treated
> as members of the SparseMatrix inheritance family.  It's bad enough
> that we can't simultaneously keep the LaspackMatrix interface
> up-to-date with changes to the Petsc one, and now we want to add one
> more family member? ;-)

Speaking of which, I forgot about the Trilinos stuff.  it seems we've
made our lives harder than necessary by having nearly all the
functions in NumericVector and SparseMatrix be pure virtual.  As long
as there is at least one pv function in the class to prevent
instantiations thereof, would it be easier on us if the rest had
default implementations of libmesh_error() and suitable a error
message?

-- 
John

On Thu, 16 Oct 2008, John Peterson wrote:

>> I think I'll withhold judgment until I've seen the new ShellMatrix
>> header file... It doesn't seem like they should necessarily be treated
>> as members of the SparseMatrix inheritance family.  It's bad enough
>> that we can't simultaneously keep the LaspackMatrix interface
>> up-to-date with changes to the Petsc one, and now we want to add one
>> more family member? ;-)

Eh; I'd rather have a nice OOP API that isn't completely implemented
than a non-OOP API that isn't completeable.

> Speaking of which, I forgot about the Trilinos stuff.  it seems we've
> made our lives harder than necessary by having nearly all the
> functions in NumericVector and SparseMatrix be pure virtual.  As long
> as there is at least one pv function in the class to prevent
> instantiations thereof, would it be easier on us if the rest had
> default implementations of libmesh_error() and suitable a error
> message?

Like libmesh_not_implemented(), you mean?  Yeah, that's probably the
best way to go.  Ideally, when you discover that your new pure virtual
method is uncompilable until you add definitions for all the
subclasses that might use it, that's just motivation to complete the
work.  In reality, none of us are that diligent, and if we're going to
be lazy anyway we ought to make laziness more convenient.  ;-)
---
Roy

> In reality, none of us are that diligent, and if we're going to
> be lazy anyway we ought to make laziness more convenient.  ;-

Agree 100%.

-Ben

Dear John/Roy/all

On Thu, 16 Oct 2008, John Peterson wrote:

> On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner <roystgnr@...> wrote:
>>
>> On Thu, 16 Oct 2008, Tim Kroeger wrote:
>>
>> I'm a bit unhappy with the current inheritance tree - there's no
>> inheritance between ShellMatrix and SparseMatrix, even though the
>> latter is already practically an abstract base class that declares the
>> functions the former will need to use?  There's no getting around the
>> need to duplicate a lot of function definitions (get_diagonal() is
>> going to be different in the different classes no matter how things
>> are factored) but it seems like a little refactoring (make both
>> ShellMatrix and SparseMatrix inherit from a MatrixBase?) would at
>> least give us more concise code using them.

That sounds right for me in a way, but on the other hand, I'm not 
quite sure.  I mean, it is difficult to decide which methods to 
include in that class and which not.  The first attempt would be to 
include exactly those that are common to SparseMatrix and ShellMatrix 
(BTW, what about DenseMatrix?), but perhaps somebody wants to 
implement another type of matrix that has again more or less common 
properties with the matrices implemented so far.  If you don't mind, I 
would leave ShellMatrix as completely disjoint from SparseMatrix for 
the moment.  Of course, you may feel free to combine them in a way 
once you included them into you repository.

> I think I'll withhold judgment until I've seen the new ShellMatrix
> header file...

Oops, it ought to be included in the patch, but I forgot that "svn 
diff" wouldn't care about it unless I did "svn add" before.  Please 
find attached now a hopefully complete patch.

(Otherwise, the patch is the same as last time because I didn't have 
time for that today.)

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Fri, Oct 17, 2008 at 7:17 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear John/Roy/all
>
> On Thu, 16 Oct 2008, John Peterson wrote:
>
>> On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner <roystgnr@...>
>> wrote:
>>>
>>> On Thu, 16 Oct 2008, Tim Kroeger wrote:
>>>
>>> I'm a bit unhappy with the current inheritance tree - there's no
>>> inheritance between ShellMatrix and SparseMatrix, even though the
>>> latter is already practically an abstract base class that declares the
>>> functions the former will need to use?  There's no getting around the
>>> need to duplicate a lot of function definitions (get_diagonal() is
>>> going to be different in the different classes no matter how things
>>> are factored) but it seems like a little refactoring (make both
>>> ShellMatrix and SparseMatrix inherit from a MatrixBase?) would at
>>> least give us more concise code using them.
>
> That sounds right for me in a way, but on the other hand, I'm not quite
> sure.  I mean, it is difficult to decide which methods to include in that
> class and which not.  The first attempt would be to include exactly those
> that are common to SparseMatrix and ShellMatrix (BTW, what about
> DenseMatrix?), but perhaps somebody wants to implement another type of
> matrix that has again more or less common properties with the matrices
> implemented so far.  If you don't mind, I would leave ShellMatrix as
> completely disjoint from SparseMatrix for the moment.  Of course, you may
> feel free to combine them in a way once you included them into you
> repository.
>
>> I think I'll withhold judgment until I've seen the new ShellMatrix
>> header file...

After looking at the ShellMatrix instantiations from the patch, I'm
now leaning towards Roy's suggestion that ShellMatrix should
eventually be derived from SparseMatrix.  I think adding an extra base
class (eg. MatrixBase) for sparse matrices is unnecessary and maybe a
bit confusing.  Making all but one or two (how about m() and n()?) of
the currently-pure-virtual functions in the SparseMatrix interface be
calls to not_implemented() will help us a lot out in this endeavour.

Anyway, my main reasons for thinking this way now are:

1.) The current ShellMatrix interface is only about five functions:
m(), n(), vector_mult(), vector_mult_add(), and get_diagonal().  But
if all goes well we may eventually want to, say, print a ShellMatrix
or get its L1-norm or do any of the umpteen other capabilities
currently supported by SparseMatrix.  So this reason is really geared
towards future-proofing the design.

2.) We would get the standard OOP capability that allows existing
library routines taking a SparseMatrix reference or pointer to take a
ShellMatrix pointer without other changes to the code!

3.) The SparseShellMatrix from Tim's patch is a great example of the
"Proxy" Design Pattern (http://en.wikipedia.org/wiki/Proxy_pattern).
In this pattern, both the proxy (SparseShellMatrix) and the the real
subject (PetscMatrix) are generally derived from the same base
(SparseMatrix) and the proxy's functionality is implemented in terms
of the real subject.  Deriving ShellMatrix (and therefore
SparseShellMatrix) from SparseMatrix makes this relationship clear.

I'm not stipulating this as a requirement for checking in the initial
patch, of course.  I volunteer to work on unifying the
SparseMatrix/ShellMatrix hierarchy around once Tim's patch has been
checked in.

-- 
John

Dear John,

Thank you for your comments.

If you like to unify SparseMatrix and ShellMatrix after commiting my 
patch, I would find that completely agreeable.  I would guess that 
it's a lot more difficult than it looks like on first glance, but I 
might be wrong.

One problem would be that this makes it more difficult to 
differentiate (e.g. in LinearSolver::solve()) which type of matrix is 
actually present.  On the other hand, that might easily be solvable 
using some dynamic_cast statements (which I am not too familiar with).

Also, I would suggest to derive SparseMatrix from ShellMatrix rather 
than vice versa.  This could actually make SparseShellMatrix 
unnecessary.  This is of course an advantage of the unification.


By the way, my application seems to get incorrect results at the 
moment, and I have to find out whether the error is on the ShellMatrix 
side or on the application side.  Since every test run takes several 
hours, you should have some more patience.  Also, the cluster is still 
not usable, so that I cannot test in parallel.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Mon, Oct 20, 2008 at 8:01 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear John,
>
> One problem would be that this makes it more difficult to differentiate
> (e.g. in LinearSolver::solve()) which type of matrix is actually present.
>  On the other hand, that might easily be solvable using some dynamic_cast
> statements (which I am not too familiar with).

Yeah, overload resolution is done at compile time and it may be
ambiguous or just choose the first solve() function which works.
We'll figure that out when we get there, if worse comes to worst we
might have to, as you say, do some dynamic cast tests within the solve
function itself, though I hope to avoid that.

Another possibility might be to subclass the PetscLinearSolver or
create a ShellLinearSolver which implements the desired functionality.
 There's no reason to completely tie the "shell" concept to the Petsc
linear solver, even though that is the only good example we currently
have of one.


> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than
> vice versa.  This could actually make SparseShellMatrix unnecessary.  This
> is of course an advantage of the unification.

This one I disagree with.  Inheritance should always follow the "is a"
organizational semantic, and a SparseMatrix (as we now define it) is
most definitely *not* a ShellMatrix.

-- 
John

On Mon, 20 Oct 2008, John Peterson wrote:

> On Mon, Oct 20, 2008 at 8:01 AM, Tim Kroeger
> <tim.kroeger@...> wrote:
>> Dear John,
>>
>> One problem would be that this makes it more difficult to differentiate
>> (e.g. in LinearSolver::solve()) which type of matrix is actually present.
>>  On the other hand, that might easily be solvable using some dynamic_cast
>> statements (which I am not too familiar with).
>
> Yeah, overload resolution is done at compile time and it may be
> ambiguous or just choose the first solve() function which works.
> We'll figure that out when we get there, if worse comes to worst we
> might have to, as you say, do some dynamic cast tests within the solve
> function itself, though I hope to avoid that.

There's no way around having to do a dynamic cast test in solve(), I
don't think.  I'd be happy to be proven wrong, though.

But anyway, in the long run it's better to have a nice clear interface
with the junk like dynamic_cast encapsulated away behind the scenes.

>> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than
>> vice versa.  This could actually make SparseShellMatrix unnecessary.  This
>> is of course an advantage of the unification.
>
> This one I disagree with.  Inheritance should always follow the "is a"
> organizational semantic, and a SparseMatrix (as we now define it) is
> most definitely *not* a ShellMatrix.

Yes, but both of them are a MatrixBase.

Yes, it is easier to add new classes when someone else has already
volunteered to do the work of writing them.  ;-)
---
Roy

On Mon, Oct 20, 2008 at 11:14 AM, Roy Stogner <roystgnr@...> wrote:
>>> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than
>>> vice versa.  This could actually make SparseShellMatrix unnecessary.
>>>  This
>>> is of course an advantage of the unification.
>>
>> This one I disagree with.  Inheritance should always follow the "is a"
>> organizational semantic, and a SparseMatrix (as we now define it) is
>> most definitely *not* a ShellMatrix.
>
> Yes, but both of them are a MatrixBase.
>
> Yes, it is easier to add new classes when someone else has already
> volunteered to do the work of writing them.  ;-)

In all seriousness, I still don't see the compelling argument for this
extra base class.  You still have to do dynamic casting in the solve
function, for example.

-- 
John

On Mon, 20 Oct 2008, John Peterson wrote:

> On Mon, Oct 20, 2008 at 11:14 AM, Roy Stogner <roystgnr@...> wrote:
>>>> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than
>>>> vice versa.  This could actually make SparseShellMatrix unnecessary.
>>>>  This
>>>> is of course an advantage of the unification.
>>>
>>> This one I disagree with.  Inheritance should always follow the "is a"
>>> organizational semantic, and a SparseMatrix (as we now define it) is
>>> most definitely *not* a ShellMatrix.
>>
>> Yes, but both of them are a MatrixBase.
>>
>> Yes, it is easier to add new classes when someone else has already
>> volunteered to do the work of writing them.  ;-)
>
> In all seriousness, I still don't see the compelling argument for this
> extra base class.

Because SparseMatrix is-not-a ShellMatrix (doesn't need to attach or
even acknowledge the existance of a user's matvec shell function), and
ShellMatrix is-not-a SparseMatrix (doesn't need sparsity pattern
details in the constructor, can't return or set individual matrix
entries), but both of them are-a MatrixBase (can perform matrix-vector
products on NumericVectors, at least, and with some tricks we can pass
them both to what are conceptually the same solvers, scale them, add
them together, etc).
---
Roy

On Oct 20, 2008, at 10:23 AM, Roy Stogner wrote:
> Because SparseMatrix is-not-a ShellMatrix (doesn't need to attach or
> even acknowledge the existance of a user's matvec shell function), and
> ShellMatrix is-not-a SparseMatrix (doesn't need sparsity pattern
> details in the constructor, can't return or set individual matrix
> entries), but both of them are-a MatrixBase (can perform matrix-vector
> products on NumericVectors, at least, and with some tricks we can pass
> them both to what are conceptually the same solvers, scale them, add
> them together, etc).

FWIW - I completely agree with this.

BTW: Is this another discussion that should get moved over to libmesh- 
devel?

Derek

Dear John,

On Mon, 20 Oct 2008, John Peterson wrote:

>> Also, I would suggest to derive SparseMatrix from ShellMatrix rather than
>> vice versa.  This could actually make SparseShellMatrix unnecessary.  This
>> is of course an advantage of the unification.
>
> This one I disagree with.  Inheritance should always follow the "is a"
> organizational semantic, and a SparseMatrix (as we now define it) is
> most definitely *not* a ShellMatrix.

Nor is a ShellMatrix a SparseMatrix.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

Dear John,

On Tue, 21 Oct 2008, Tim Kroeger wrote:

> On Mon, 20 Oct 2008, John Peterson wrote:
>
>> This one I disagree with.  Inheritance should always follow the "is a"
>> organizational semantic, and a SparseMatrix (as we now define it) is
>> most definitely *not* a ShellMatrix.
>
> Nor is a ShellMatrix a SparseMatrix.

What I mean is that SparseMatrix has all the functionality that 
ShellMatrix requires plus a lot of extra functionality.  If we had a 
MatrixBase as Roy suggests, then what would be the additional 
functionality of ShellMatrix with respect to MatrixBase?  I can't see 
any.  Hence, what about renaming ShellMatrix to MatrixBase and then 
deriving SparseMatrix from that class?

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

Dear John/Roy/all,

Another question pops up in connection with my ShellMatrix, in 
particular with the vector v and w that form a TensorShellMatrix. 
When I assemble v and w (where the matrix B = v w'), I certainly have 
to call something similar to DofMap::constrain_element_matrix(), but 
I'm not quite sure what the correct thing is, since I'm not familiar 
with libMesh's constaining mechanism.

My guess would be to call

DofMap::constrain_element_vector(v,dof_indices);

and do nothing to w.  Can you please let me know whether this is 
correct?

Also, I understand that all of these constraining functions 
occassionally have to change their dof_indices argument, hence if I 
assemble a SparseMatrix and a TensorShellMatrix at the same time, the 
following is obviously wrong:

DofMap::constrain_element_matrix_and_vector(Ke,Fe,dof_indices);
DofMap::constrain_element_vector(v,dof_indices);

I mean, the second call will get the modified dof_indices and 
interpret v in the wrong way.  Hence, I guess the correct thing to do 
is to implement a new method 
DofMap::constrain_element_matrix_and_two_vectors() -- or even 
DofMap::constrain_element_matrix_and_several_vectors() (taking a 
vector of vectors).

Please let me know what you suggest.  I think I would be able to 
implement such methods myself.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen

On Tue, Oct 21, 2008 at 1:54 AM, Tim Kroeger
<tim.kroeger@...> wrote:
> Dear John,
>
> What I mean is that SparseMatrix has all the functionality that ShellMatrix
> requires plus a lot of extra functionality.  If we had a MatrixBase as Roy
> suggests, then what would be the additional functionality of ShellMatrix
> with respect to MatrixBase?  I can't see any.  Hence, what about renaming
> ShellMatrix to MatrixBase and then deriving SparseMatrix from that class?

OK.  I'd also like to change the name TensorShellMatrix ->
DyadShellMatrix (http://en.wikipedia.org/wiki/Dyadic_product), if
that's alright.  The attached updated class tree shows how I now
envision the hierarchy going.

I decided to employ diamond inheritance for the PetscMatrix because it
was the neatest way to resolve the simultaneous need to have a single
base for all Petsc-type matrices (so that PetscLinearSolver routines
can work seamlessly) while still allowing shell matrices to come from
a separate family tree than SparseMatrices come from.  We'll just be
careful to use the proper virtual inheritance.

The MatCreateShell and MatShellSetOperation calls in Tim's last patch
will move from the PetscLinearSolver class to the new PetscShellMatrix
class, which shall be derived from PetscMatrixBase.  PetscLinearSolver
code will then dynamic_cast internally to PetscMatrixBase* to get at
the 'Mat' object.  Generic LinearSolver routines can take
NumericMatrix<T>& references instead of SparseMatrix<T>& references.
I'd be curious to hear about other ideas/suggestions...

-- 
John

Dear John,

On Tue, 21 Oct 2008, John Peterson wrote:

> On Tue, Oct 21, 2008 at 1:54 AM, Tim Kroeger <tim.kroeger@...> wrote:
>>
>> What I mean is that SparseMatrix has all the functionality that ShellMatrix
>> requires plus a lot of extra functionality.  If we had a MatrixBase as Roy
>> suggests, then what would be the additional functionality of ShellMatrix
>> with respect to MatrixBase?  I can't see any.  Hence, what about renaming
>> ShellMatrix to MatrixBase and then deriving SparseMatrix from that class?
>
> OK.  I'd also like to change the name TensorShellMatrix ->
> DyadShellMatrix (http://en.wikipedia.org/wiki/Dyadic_product), if
> that's alright.

Yes, of course.

> The attached updated class tree shows how I now
> envision the hierarchy going.

Your inheritance hierarchy makes user code that uses the shell 
matrices functionality solver dependent: The user has to explicitly 
instantiate e.g. PetscDyadMatrix.  You might later implement shell 
matrices for other solver packages, and I would find it more 
appropriate to have a DyadMatrix class as the user API and let 
everything that depends on the solver package be hidden from the user.

Also, I'm somehow skeptic about the diamond inheritance, but that 
might just be because I'm not experienced in this.  Actually, 
resolving the problem that I stated above might require more diamond 
inheritance if you stick to the concept that PETSc's Mat object should 
be managed by the matrix class rather than the solver class.

Anyway, I should emphasize once more that you may of course feel free 
to restructure everything as you like once it's in the repository, 
provided that it doesn't break my application.  (API changes are of 
course also okay at this stage if I get informed what I have to 
change in my application.)

By the way: What about my question concerning matrix/vector 
contraining?  I didn't receive any answer to that from any of you 
guys.

By the way: I attached a new version of my patch.  I have added the 
possibility to use a SparseMatrix as the preconditioning matrix when 
the system matrix is a ShellMatrix.

Best Regards,

Tim

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

MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany

Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen