[PETSC #13812] Feature request: "long double" support

[Barry Smith]

  Roy,

Number of Newton iterations = 4
lid velocity = 0.0204082, prandtl # = 1, grashof # = 1
  0 SNES Function norm 1.461936698223e-01 
  1 SNES Function norm 2.507018921541e-05 
  2 SNES Function norm 5.123354965471e-13 
  3 SNES Function norm 1.676580470503e-28 
  4 SNES Function norm 5.094694802293e-32 
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE

I have added preliminary support for compiling PETSc using
long doubles (128 bit on some systems like the Apple). The 
config/configure.py flag is --with-precision=longdouble

Current limitations:
  1) have not added support for using Fortran (I think this
     is doable assuming they use the same floating point
     representation)
  2) no support for binary IO (this is doable, just need
     work arrays, etc. like the code for --with-precision=float)
  3) must use --download-c-blas-lapack (assuming you have native
     BLAS/LAPACK for 128 bit (I think IBM does) we could add
     support to use it.)
  4) some systems don't support %Lg for ASCII formating (for example,
     Apple doesn't support it, but it does support long double :-().
  5) no support for complex

  I won't work on 1,2 or 3 until there is a user who needs it.

   Barry

On Fri, 14 Oct 2005, Roy Stogner wrote:

> On Fri, 14 Oct 2005, Barry Smith wrote:
> 
> >   Yes I think it would be reasonable, and maybe possible :-).
> > 
> >   The kicker will be the BLAS calls, we'd either need to 1) write the
> > name translations if your system has some long double BLAS
> 
> I'm afraid it doesn't.  I'm using the BLAS provided with Fedora Core,
> which has s/d/c/z versions of every call and no more.  I've seen
> proposals to add versions "k" for long double and "q" for long double
> complex, but I don't know if this has been coded or widely adopted
> yet.
> 
> > or 2) compile f2cblaslapack or blaslapack in long double or
> 
> I'd be worried about portability between F77 and C code.  I know on
> gcc/x86 for example the ABI for 80 bit long doubles can be padded out
> to 96 bits or 128 bits based on compiler flags.  I'm not sure if the
> 96 bit version corresponds to Fortran's REAL*12, and I don't think the
> 128 bit version corresponds to REAL*16.
> 
> You could probably start with the double f2cblaslapack routines and
> get working code out of 's/double/long double/g', though.  Trying to
> maintain code that came out of a wizard or autotranslator is a
> nightmare, but it should be fine for code that doesn't require a lot
> of updates.  This is probably the way to go.
> 
> > 3) turn on PETSc with Fortran kernels to replace all the BLAS (and
> > compile the Fortran with long double).
> 
> This wouldn't work unless I'm wrong about the incompatibility between
> Fortran and C/C++ extended precision datatypes.
> 
> > None of these are theoretically difficult, but the process could be
> > a little messy.
> 
> That's about what I thought.  If it wasn't going to be messy I would
> have sent you a patch, not a feature request.  ;-)
> 
> >   Changes to the rest of PETSc would be (in theory) trivial.
> 
> That's what I hoped, but there's so many places in PETSc where
> "double" is used literally rather than indirectly through PetscScalar;
> I don't know just how much of that would have to be extended.
> 
> For a code you're already familiar with it's not hard.  I've been
> checking libMesh for float and long double compatibility, and the only
> unexpected problems were some *printf format strings and tolerance
> constants that need to be precision dependent.  The only hard problem
> is interfacing with outside libraries - it might be reasonable to
> throw away precision and truncate to double before interacing with I/O
> code, but if there's a third party library handling intermediate
> calculations where you need that precision, then the only solution is,
> well, to pester their email address like this.
> 
> Thanks,
> ---
> Roy Stogner
> 
>