The biggest reason we haven't added "multi-vectors" is more design decisions that have to be made then actual code implementation. There are two important design decisions: I) Introduce a new class or overload Vec? 1) Introduce a Vecs class, essentially copying the interface from Vec and implement the methods. A Vecs is essentially a collection of Vec. Then provide methods such as MatMults(), MatSolves(), PCApplys(), PCSetUps(), KSPSolves() etc 2) Overload the current Vec class so that a Vec is either a "regular" Vec or a collection of vectors. Modify the Vec interface and methods as needed, then modify the MatMult() etc as needed 3) Do 1) and have it all working well and then use this to quickly do 2) and throw away the seperate concept of Vecs (without losing the fastest possible performance for a "new" Vec that has a single vector). The advantage of 1) is we never break the current interface. 2) has the advantage of not introducing a new class; and since software comprehension difficulty grows exponentially with the number of classes this is a good thing. The advantage of 3) is we get an incremental approach to 2) (incremental is often good). Disadvantage of 3) is that users have to live with first the Vecs interface (but on the other hand most people won't use Vecs anyways) and then go back to Vec. I think 3) is the correct approach. II) Default storage of "multivectors"? (note: the interfaces for Vecs will, of course, be data storage format neutral, but the code we actually write cannot be) 1) Have the Vecs be essentially an array of Vec. 2) Have a Vecs essentially be a double array where the "values" of the "vectors" are interlaced. The advantage of 1) is that it is easier to implement and doing things like deflation (removing one of the vectors in a Vec) doesn't require any copies of the vector data. The advantage of 2) is better performance, and since we are using block for better performance we should take advantage of the best possible performance. I think implementing both 1) and 2) is appropriate, starting with 1) since getting it completed would be faster. Barry Note: the Vecs currently in petscvec.h is just a toy and can be ignored. On Fri, 7 Apr 2006, Victor Eijkhout wrote:
The question has come up before. If all goes according to plan, then Bill Gropp will hire a student for the summer to implement it. On a project that I lead, and that has use for multivectors for a different purpose, there is actually budget. Feel free to recommend someone to Bill if you know of a good candidate. Does your former advisor have students who are familiar with the matter?
Victor.
On Apr 7, 2006, at 1:05 PM, Richard Tran Mills wrote:
Hi Guys,
My former Ph.D. advisor, Andreas Stathopoulos ([email protected]), has released an eigensolver package, PRIMME, that he is interested in interfacing as an external package with PETSc. I told him a bit about how he could do this, but what I couldn't really answer was his questions about "multivector" support in PETSc. PRIMME implements block Davidson/Jacobi-Davidson type methods, and thus needs to do a lot of operations like sparse matrix-multivector multiplies. I haven't seen anything like this in PETSc (and since PETSc doesn't have block Krylov solvers, I suppose there is currently no need). Has anyone done any work on stuff like this, or does anyone have any plans to? Any thoughts on this?
Cheers, Richard