<div>I saw a good talk by Chris Baker at SIAM about multiprecision solvers, and was suitably embarrassed that we cannot do</div><div>this. I still strongly believe that templates are the wrong solution, and mucking up perfectly good interfaces with a</div>
<div>bunch of types we do not care about is wrong. I now have a proposal, outlined below.</div><div><br></div><div>We start with the idea of a solver with a certain precision, so that each solver can report the precision it computes</div>
<div>in. Let's restrict ourselves to solvers that only do simple vector operations for the time being. For this then, we</div><div>need a Vec object which stores values and computes with this precision, which will be discussed below. Upon entering,</div>
<div>we call VecConvert() on the input vector to get the correct precision, which just references the object if no conversion</div><div>is necessary. We VecConvert() the output vector as well before filling it, and the VecConvert() to the output on the way</div>
<div>out. If not conversion is necessary, we get the ssme pointer for the first convert, and no-op for the second. All</div><div>internal vectors are created with the solver precision. This all depends on having Vec objects of different precisions.</div>
<div><br></div><div>I propose we have a different Vec implementation for each precision, maybe saving code by having configure or make</div><div>generate the different precisions from a single template. Most Vec operations work by first calling VecGetArray() and</div>
<div>then some BLAS calls. I would now have VecGetArray() do copy-in/copy-out if internal storage has a different precision</div><div>than PetscScalar, and add the function VecGetArrayPrecise(Vec v, PetscDataType precision, void **array) which returns</div>
<div>the internal storage if it is the correct precision, and fails otherwise. If we switch implementations to</div><div>VecGetArrayPrecise(), then vectors will work with other vectors of the same precision and fail otherwise. I think this</div>
<div>is the correct level of sanity with multi-precision code (demand explicit conversions of dissimilar types).</div><div><br></div><div>I am only really interested in trying this with float+double (Chris spent all of his time on qd128 which I think has</div>
<div>no real scientific relevance, but of course we will do this too). The code changes are relatively minor:</div><div>  - Add a precision (PetscDataType) to Vec</div><div>  - One new Vec method, VecGetArrayPrecise()</div>
<div>  - Make VecDuplicate() respect precision</div><div>  - New Vec implementation for float, built from a template by Python</div><div>If we give a float Vec to Picard, it should compute everything in float since the Vecs are duplicated. We can add the</div>
<div>solver convert code later if everything works as expected. I think the first interesting test would be a nonlinear float</div><div>preconditioner for a double SNES.</div><div><br></div><div>I am not interested in being able to use whatever type I want everywhere. That will never work, and its not clear there is</div>
<div>any benefit anyway. What I want is just to have solvers of different precisions work together with explicit conversion at</div><div>the boundary.</div><div><br></div><div>Why will this not work? If it will, why will it not be as easy as stated above?</div>
<div><br></div><div>    Matt</div><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>