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
