[petsc-dev] Multiprecision solver composition

[Matthew Knepley]

On Sat, Jul 14, 2012 at 2:57 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>    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.
>

I am fine with all the regular functions to take PetscScalar. That keeps
the interface simple, and I don't think you can make
a single argument against casting around single scalars or small arrays
like VecMDot().


>    VecGetArray() requesting the precision type may be useful as a starting
> point.


I think the only ones we need to be careful with are the direct value
access, MatSetValues(), VecGetArray(), etc. For this, I
recommend what I said, have another function which takes an explicit
precision, and have the normal method do automatic
conversion.

If there are no objections more damaging than that, I think we might as
well try it.

   Matt


>
>    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
>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20120714/a60ccde4/attachment.html>