It's already there! :-) From petscsys.h
/*EC
PetscPrecision - indicates what precision the object is using. This is
currently not used.
Level: advanced
.seealso: PetscObjectSetPrecision()
E*/
typedef enum { PETSC_PRECISION_SINGLE=4,PETSC_PRECISION_DOUBLE=8 }
PetscPrecision;
PETSC_EXTERN const char *PetscPrecisions[];
Matt,
The biggest stumbling block I've had with doing this in C is not the part
you talk about below but the access like VecSetValues(), VecAXPY() etc (and
all the Mat ones) that
take a scale argument or worse array of scales. How to manage that? You can
pass any precision into any object, or only the precision of the object? Always
double precision for these arguments? but that is too limiting .... It is hard
for me to see how to prevent an explosion of "templated" interfaces of some
kind or another.
VecGetArray() requesting the precision type may be useful as a starting
point.
Barry
On Jul 14, 2012, at 10:09 AM, Matthew Knepley wrote:
> I saw a good talk by Chris Baker at SIAM about multiprecision solvers, and
> was suitably embarrassed that we cannot do
> this. I still strongly believe that templates are the wrong solution, and
> mucking up perfectly good interfaces with a
> bunch of types we do not care about is wrong. I now have a proposal, outlined
> below.
>
> We start with the idea of a solver with a certain precision, so that each
> solver can report the precision it computes
> in. Let's restrict ourselves to solvers that only do simple vector operations
> for the time being. For this then, we
> need a Vec object which stores values and computes with this precision, which
> will be discussed below. Upon entering,
> we call VecConvert() on the input vector to get the correct precision, which
> just references the object if no conversion
> is necessary. We VecConvert() the output vector as well before filling it,
> and the VecConvert() to the output on the way
> out. If not conversion is necessary, we get the ssme pointer for the first
> convert, and no-op for the second. All
> internal vectors are created with the solver precision. This all depends on
> having Vec objects of different precisions.
>
> I propose we have a different Vec implementation for each precision, maybe
> saving code by having configure or make
> generate the different precisions from a single template. Most Vec operations
> work by first calling VecGetArray() and
> then some BLAS calls. I would now have VecGetArray() do copy-in/copy-out if
> internal storage has a different precision
> than PetscScalar, and add the function VecGetArrayPrecise(Vec v,
> PetscDataType precision, void **array) which returns
> the internal storage if it is the correct precision, and fails otherwise. If
> we switch implementations to
> VecGetArrayPrecise(), then vectors will work with other vectors of the same
> precision and fail otherwise. I think this
> is the correct level of sanity with multi-precision code (demand explicit
> conversions of dissimilar types).
>
> I am only really interested in trying this with float+double (Chris spent all
> of his time on qd128 which I think has
> no real scientific relevance, but of course we will do this too). The code
> changes are relatively minor:
> - Add a precision (PetscDataType) to Vec
> - One new Vec method, VecGetArrayPrecise()
> - Make VecDuplicate() respect precision
> - New Vec implementation for float, built from a template by Python
> If we give a float Vec to Picard, it should compute everything in float since
> the Vecs are duplicated. We can add the
> solver convert code later if everything works as expected. I think the first
> interesting test would be a nonlinear float
> preconditioner for a double SNES.
>
> I am not interested in being able to use whatever type I want everywhere.
> That will never work, and its not clear there is
> any benefit anyway. What I want is just to have solvers of different
> precisions work together with explicit conversion at
> the boundary.
>
> Why will this not work? If it will, why will it not be as easy as stated
> above?
>
> Matt
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
