On Tue, 7 Oct 2008, John Peterson wrote: > On Tue, Oct 7, 2008 at 10:32 AM, Tim Kroeger > <[EMAIL PROTECTED]> 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 ------------------------------------------------------------------------- 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 [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-users
