[petsc-dev] Simultaneous use of different precisions

[Barry Smith]

  Jed,

    It is very experimental, it may be a complete disaster and I may need to remove some of the added code so don't go assuming there is any clear sailing about it.

    I've been trying to figure out how to "support different precisions" (where even that term support is ill-defined) for over ten years within the context of PETSc objects and still don't have a handle on how to do it. 

    Reason to want the support: (1) it is crazy to always solve our linear systems in double precision in the context of nonlinear solves - wastes lots of memory and lots of time. (2) less important to me, it would be nice to have complex numbers and real numbers in the same code.

    Well, why not just use C++ templates to support it (like everyone else)? (1) I hate beyond words the standard C++ template model where one would have Vec<double> etc to indicate that the Vec is templated on double. I want to write Vec; the data is encapsulated inside it, it completely breaks the encapsulation to expose the storage type in the declaration. (2) As soon as anyone starts using C++ templates they sink inside this hole of utter incomprehensibility where no one understands the codes and compilers generate worthless error messages. I submit almost no users would tolerate a gaudy templated up PETSc, not even me because it throws out data enscapsulation, plus throws out Fortran. Now, I could be totally wrong about this, wish I were and that it was possible to do a "simple" templated PETSc with everything we want, please, please, please show me how to do it.

  What do I mean by "support different precisions"? Well ideally (call this Model I) it would mean I can declare an object to use a precision internally (like double) and then use it anywhere with any other object. So, for example, if I declared a Vec to be double internally I could call VecSetValues() with either float or double numbers, for example, and find its inner product with another Vec that happens to be internally either double or float. Given the current PETSc object model and source code I consider this "too hard".  Frankly I think this is too hard with any language we have today.

So instead (call this Model II)  what about being able to create double objects that can interact only with other double objects and having all calling sequences with double values (and float objects only interact with other float objects), these could be KSP, PC, Mat, Vec. In addition have a way of "copying" single Vecs to double Vecs and vs versa. (not bothering to have copies between single and double matrices). Then SNES could (in theory anyways) trivially switch between single and double linear solvers.  This, I think, can be done using ONLY simple C++ function polymorphism (allowing VecSetValues(double*) and VecSetValues(single *) for example) and since Fortran 90 can (in theory anyways) support this same simple polymorphism (in theory anyways) it may be possible to have the same user model in Fortran with appropriate wrappers (and python too I assume). 

Now I think both Model I and Model II can be handled "mostly" just by suitable code organization but Model I is much much complicated and (I think) requires mimicing templates alot (for example auto-generating all the possible variants of dot product between float and double). I am seeing how far I can get with Model II with just minimal code reorganization. Most of the organization is good anyways and won't need to be undone when I realize this is a trainwreck.

So what I need is a PetscObjectSetPrecision(), reorganized functions in source files so that single and double libraries can coexist and some way to have double or single XXXCreate_YYY() called based on their precision. This is what I am doing now.

   Barry


    
On Nov 17, 2010, at 3:09 PM, Jed Brown wrote:

> Barry, I see some commits in this direction.  Would you mind giving us a very brief summary of where you are headed with this?
> 
> Jed