Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-21 06:54 ```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 ```

[Libmesh-users] Matrix free solving Tim Kroeger <tim.kroeger@ce...>
 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-07 15:42 ```On Tue, Oct 7, 2008 at 10: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? 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-07 16:07 ```On Tue, 7 Oct 2008, John Peterson wrote: > On Tue, Oct 7, 2008 at 10:32 AM, 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 ```

 Re: [Libmesh-users] Matrix free solving From: Derek Gaston - 2008-10-07 16:23 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-08 15:52 Attachments: patch ```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* 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-08 17:01 ```On Wed, Oct 8, 2008 at 10:52 AM, 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 for inner_ctx, a pair 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-08 17:05 ```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 for inner_ctx, a > pair 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-08 17:16 ```On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner 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 for inner_ctx, a >> pair 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-08 17:29 ```On Wed, 8 Oct 2008, John Peterson wrote: > On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-08 18:55 ```On Wed, Oct 8, 2008 at 12:55 PM, Jed Brown wrote: > On Wed 2008-10-08 12:42, John Peterson wrote: >> Thanks Jed... >> >> On Wed, Oct 8, 2008 at 12:24 PM, Jed Brown 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(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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-08 21:18 ```On Wed, Oct 8, 2008 at 12:05 PM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-09 12:31 Attachments: patch ```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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-09 14:48 ```On Thu, Oct 9, 2008 at 7:04 AM, 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-10 16:28 ```John Peterson wrote: > On Thu, Oct 9, 2008 at 7:04 AM, 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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-10 16:41 Attachments: patch ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-16 15:20 Attachments: patch ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-16 15:47 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-16 15:57 ```On Thu, Oct 16, 2008 at 10:20 AM, 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(static_cast(&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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-16 16:21 ```On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-16 16:44 ```On Thu, Oct 16, 2008 at 11:20 AM, John Peterson wrote: > On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-16 16:58 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Benjamin Kirk - 2008-10-16 17:21 ```> 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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-17 12:17 Attachments: patch ```Dear John/Roy/all On Thu, 16 Oct 2008, John Peterson wrote: > On Thu, Oct 16, 2008 at 10:46 AM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-17 15:51 ```On Fri, Oct 17, 2008 at 7:17 AM, 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 >> 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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-20 13:02 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-20 16:07 ```On Mon, Oct 20, 2008 at 8:01 AM, 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-20 16:14 ```On Mon, 20 Oct 2008, John Peterson wrote: > On Mon, Oct 20, 2008 at 8:01 AM, 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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-20 16:17 ```On Mon, Oct 20, 2008 at 11:14 AM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: Roy Stogner - 2008-10-20 16:23 ```On Mon, 20 Oct 2008, John Peterson wrote: > On Mon, Oct 20, 2008 at 11:14 AM, Roy Stogner 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 ```

 Re: [Libmesh-users] Matrix free solving From: Derek Gaston - 2008-10-20 17:03 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-21 06:29 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-21 06:54 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-21 14:51 ```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 ```

 Re: [Libmesh-users] Matrix free solving From: John Peterson - 2008-10-21 22:18 ```On Tue, Oct 21, 2008 at 1:54 AM, 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& references instead of SparseMatrix& references. I'd be curious to hear about other ideas/suggestions... -- John ```

 Re: [Libmesh-users] Matrix free solving From: Tim Kroeger - 2008-10-22 10:38 Attachments: patch ```Dear John, On Tue, 21 Oct 2008, John Peterson wrote: > On Tue, Oct 21, 2008 at 1:54 AM, 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 ```