[petsc-dev] Multiprecision solver composition

[Matthew Knepley]

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