Re: [Boost-users] converting to float128 from double

Subject: Re: [Boost-users] converting to float128 from double
From: Maximilian Winkler (maximilian.winkler_at_[hidden])
Date: 2016-05-04 13:55:29

Hello,

changing line 10 (the 'open' statement) to

    open(unit=68,file='coup_VQZ+5jp_iv1.dat',status='old',action='read')

the program compiles fine for me using gfortran4.9.2.
An input file 'coup_VQZ+5jp_iv1.dat' of the form

 1       1
 1       1
 1       1
...     ...

will allow the program to complete. For 15'000 e/g pairs the output will
be
                  ...                                            ...                              ...

 1.00000000000000000000000000000000000         1.00000000000000000000000000000000000             14999
 1.00000000000000000000000000000000000         1.00000000000000000000000000000000000             15000
  
   15000.0000000000000000000000000000000      
 End of program
 Number of energy points =        15000

Best Regards,

Max

On Wed, May 04, 2016 at 04:29:11PM +0000, Cooper, Bridgette R D wrote:
> 
>         program run_stieltjes
>         implicit none
>         integer  num,printflag,unitout,io
>         integer np,i
>         parameter (np=15000)
>         real*8 e8(np),g8(np),gamma0
> 
>         io=0
>         num=0
>         open(unit=68,file='coup_VQZ+5jp_iv1.dat',status='old', readonly)
>         do i=1,15000
>                 read(68,*,iostat=io) e8(i),g8(i)
>                 if (io > 0) then
>                         write(6,*)'Check input file of correct format'
>                         exit
>                 else if (io < 0) then
>                         write(6,*)'Reach the end of file'
>                         write(6,*)'Number of energy points = ',num
>                         exit
>                 else
>                 num=num+1
>                 end if
>         end do
>         close(68)
>         open(unit=65,file='Stieltjes_VQZ+5jp_iv2.dat')
>         unitout=65
>         printflag=1
>         call stieltjes(num,e8,g8,gamma0,65,printflag)
>         write(6,*)'End of program'
>         write(6,*)'Number of energy points = ',num
>         close(65)
>         end
> 
>         subroutine stieltjes (num,e8,g8,gamma0,unitout,printflag)
>         implicit none
>         integer :: num,nmax,np,npol,unitout,printflag,i
>         parameter (np=15000,nmax=1000)
>         real*16 :: e_point(np),g_point(np)
>         real*16 :: bcoef(0:nmax)
>         real*8  :: e8(np),g8(np),gamma0
> 
>         do i=1,num
>            e_point(i)=e8(i)
>            g_point(i)=g8(i)
>            if (printflag.ne.0) write(6,*) e_point(i),g_point(i),i
>         end do
>           write(6,*)' '
>         bcoef(0)=0.q0
>         do i=1,num
>            bcoef(0)=bcoef(0)+g_point(i)
>         end do
>         write(6,*) bcoef(0)
>         return
>         end
> 
> This compiles with the intel compiler, but not gfortran unfortunately.
> ________________________________________
> From: Boost-users <boost-users-bounces_at_[hidden]> on behalf of John Maddock <jz.maddock_at_[hidden]>
> Sent: Wednesday, May 4, 2016 4:49:07 PM
> To: boost-users_at_[hidden]
> Subject: Re: [Boost-users] converting to float128 from double
> 
> On 04/05/2016 13:47, Cooper, Bridgette R D wrote:
> > Hi John,
> >
> > I tried again. I compile pretty much exactly the function you wrote:
> 
> Bridgette, without the input values I still can't go about reproducing
> here - I need a *complete* test case plus the matching FORTRAN code.
> 
> John.
> 
> >
> > #include <boost/multiprecision/float128.hpp>
> > #include <boost/math/cstdfloat/cstdfloat_types.hpp>
> > #include <boost/math/cstdfloat/cstdfloat_limits.hpp>
> > #include <boost/math/cstdfloat/cstdfloat_cmath.hpp>
> > #include <boost/math/cstdfloat/cstdfloat_iostream.hpp>
> > extern "C"{
> > #include <quadmath.h>
> > }
> > #include <iostream>
> > #include <fstream>
> > #include <sstream>
> > #include <vector>
> > #include <iterator>
> >
> >          using namespace std;
> >          using namespace boost::multiprecision;
> > double stieltjes_gamma(const std::vector<double> & e8, const std::vector<double> & g8, int printlevel)
> > {
> >      assert(e8.size() == g8.size());  // precondition?
> >      std::vector<float128> epoint(e8.begin(), e8.end()),
> > gpoint(g8.begin(), g8.end()), bcoef(e8.size());
> >
> >      bcoef[0] = 0.0Q;
> >      for(int i = 0; i < e8.size(); i++)
> >      {
> >         bcoef[0] += gpoint[i];
> >      }
> >      std::cout <<
> > std::setprecision(std::numeric_limits<float128>::max_digits10) <<
> > bcoef.at(0) << std::endl;
> >      return (0.0);
> > }
> >
> > The answer is still:
> >
> > 0.0345217724606016853085032291117978629
> >
> > Whereas the fortran gives the answer to be:
> >
> > 3.452175962187407162907246287483402E-0002
> >
> > I really appreciate you taking the time to respond to my emails. I really don't understand the discrepancy.
> >
> > I have also tried things like:
> >
> > epoint.at(i)=static_cast<float128>(e8.at(i));
> >
> > But still cannot reproduce the fortran value.
> >
> > I also have code that can reproduce the quad precision sqrt to match fortran code. Firstly, I have a c++ routine that calculates the sqrt (a^2+b^2) in quad precision without using the library. This reproduces the quad value from fortran, and the quad value of sqrt function using boost library. However, if the explicit type of a or b is not appended with Q I don't reproduce the correct answer, which I can understand. But what if I wanted to start with a double value of a and b, cast them up to new variables c d in float128, and want the square root of the quad numbers. How do I do this casting so that I don't get just the result casted up to a quad precision number, the sqrtq function is called?
> >
> > Thanks again for all your help.
> > ________________________________________
> > From: Boost-users <boost-users-bounces_at_[hidden]> on behalf of John Maddock <jz.maddock_at_[hidden]>
> > Sent: Wednesday, May 4, 2016 12:29:31 PM
> > To: boost-users_at_[hidden]
> > Subject: Re: [Boost-users] converting to float128 from double
> >
> > On 03/05/2016 20:58, Cooper, Bridgette R D wrote:
> >> Hi John,
> >>
> >> No, my code is several hundreds of lines long. I extracted what I thought would be quite self contained, but yeah you are right about the return type mis-match.
> >>
> >> I am not sure that this helps. What I really want is to do a lot of computation in float128 precision, and for some reason when I cast up a double vector to a float128 vector  all the operations on the float128 vector are happening a double precision, and the cast up only happens at the final stage.
> >>
> >> I've changed the initiation of epoint and gpoint vectors to the same as your code snippet and still the same problem.
> >>
> >> So for example, if I accumulate the original double vector without casting the values to float128 type and accumulate into a float128 variable, the cast up only happens to the result. This is the same number that I see as the result of accumulating the float128vector values.
> >>
> >> Hope this clarifies things a bit
> > Not really :(
> >
> > I think you're going to have to try to produce a self contained test
> > case - my guess is you'll spot the error in the process of doing that?
> >
> > Best, John.
> >
> >> ________________________________________
> >> From: Boost-users <boost-users-bounces_at_[hidden]> on behalf of John Maddock <jz.maddock_at_[hidden]>
> >> Sent: Tuesday, May 3, 2016 6:27:39 PM
> >> To: boost-users_at_[hidden]
> >> Subject: Re: [Boost-users] converting to float128 from double
> >>
> >> On 03/05/2016 17:35, Cooper, Bridgette R D wrote:
> >>> Hi John,
> >>>
> >>> Thanks for the response.
> >>>
> >>> Here are some more snippets from my code.
> >> I don't see any obvious errors, but in any case you must have extracted
> >> that from something else but it doesn't compile!
> >>
> >> Not least you can't implicitly return from a float128 to a double in
> >>
> >>        return(bcoef.at(0))
> >>
> >> So if that's compiling for you, then I suspect your declaration of bcoef.
> >>
> >> If I re-write and simplify in a slightly more C++ idiom how does it
> >> compare to your fortran code now:
> >>
> >> double stieltjes_gamma(const std::vector<double> & e8, const
> >> std::vector<double> & g8, int printlevel)
> >> {
> >>       assert(e8.size() == g8.size());  // precondition?
> >>       std::vector<float128> epoint(e8.begin(), e8.end()),
> >> gpoint(g8.begin(), g8.end()), bcoef(e8.size());
> >>
> >>       bcoef[0] = 0.0Q;
> >>       for(int i = 0; i < e8.size(); i++)
> >>       {
> >>          bcoef[0] += gpoint[i];
> >>       }
> >>       std::cout <<
> >> std::setprecision(std::numeric_limits<float128>::max_digits10) <<
> >> bcoef.at(0) << std::endl;
> >>       return static_cast<double>(bcoef[0]);
> >> }
> >>
> >> And finally... note that summing at higher precision can not actually
> >> protect you from cancellation error (in any language) since the error is
> >> inherent in the input values.
> >>
> >> HTH, John.
> >>
> >>
> >>> #include <boost/multiprecision/float128.hpp>
> >>> #include <boost/math/cstdfloat/cstdfloat_types.hpp>
> >>> #include <boost/math/cstdfloat/cstdfloat_limits.hpp>
> >>> #include <boost/math/cstdfloat/cstdfloat_cmath.hpp>
> >>> #include <boost/math/cstdfloat/cstdfloat_iostream.hpp>
> >>> #include <boost/multiprecision/detail/float_string_cvt.hpp>
> >>> #include <boost/multiprecision/detail/generic_interconvert.hpp>
> >>> extern "C"{
> >>> #include <quadmath.h>
> >>> }
> >>> #include <iostream>
> >>> #include <sstream>
> >>> #include <vector>
> >>> #include <iterator>
> >>>
> >>>            using namespace boost::multi precision;
> >>>
> >>>
> >>>     double stieltjes_gamma(int num, const std::vector<double> & e8, const std::vector<double> & g8, int printlevel )
> >>> {
> >>>
> >>>           std::vector<float128> epoint(num), gpoint(num), bcoef(num);
> >>>
> >>>
> >>>
> >>>           for (int i=0; i < num; i++)
> >>>                    {
> >>>                    epoint.at(i)=(float128(e8.at(i)));
> >>>                    gpoint.at(i)=(float128(g8.at(i)));
> >>>                    }
> >>>
> >>>     bcoef.at(0)=0.0Q;
> >>>            for (int i = 0; i < num; i++)
> >>>                    {
> >>>                    bcoef.at(0)+=gpoint.at(i)
> >>>                    }
> >>>           std::cout << std::setprecision(std::numeric_limits<float128>::max_digits10)  << bcoef.at(0)<< std::endl;
> >>>        return(bcoef.at(0))
> >>> }
> >>>
> >>> When I compare this to fortran code that does the same I don't see the same value for bcoef.at(0).
> >>>
> >>> 0.0345217724606016853085032291117978629  <- from above c++ code
> >>> 3.452175962187407162907246287483402E-0002   <- Fortran correctly utilising quad precision on same input vector
> >>>
> >>> The input vectors in this case are 7976 so errors in precision are really quickly apparrant. If I print  the gpoint vector, it gets filled with the same garbage for additional precision from double  as is for the fortran.
> >>>
> >>> Thanks,
> >>> ________________________________________
> >>> From: Boost-users <boost-users-bounces_at_[hidden]> on behalf of John Maddock <jz.maddock_at_[hidden]>
> >>> Sent: Tuesday, May 3, 2016 5:06:47 PM
> >>> To: boost-users_at_[hidden]
> >>> Subject: Re: [Boost-users] converting to float128 from double
> >>>
> >>> On 03/05/2016 13:44, Cooper, Bridgette R D wrote:
> >>>> Hi,
> >>>>
> >>>>
> >>>> I have a function that takes as an argument a vector of doubles and
> >>>> tries to convert it to a vector of float128.
> >>>>
> >>>>
> >>>> When I accumulate the values in the float128 vector it looks like it
> >>>> does double additions and only casts up at the last step.
> >>>>
> >>>>
> >>>> The initial vector is defined as conststd::vector<double>
> >>>>
> >>>> The float128 vector is defined as std::vector<float128> epoint(num)
> >>>>
> >>>>
> >>>> And I try to do the casting (inside a loop) via:
> >>>> epoint.at(i)=(float128(e8.at(i)));
> >>>>
> >>>>
> >>>> If I do the accumulation on the original vector with no casting into a
> >>>> float128, I get the same as accumulating the vector that should be of
> >>>> float128 type.
> >>>>
> >>>>
> >>>> How should I be doing the casting?
> >>>>
> >>> Sorry but you'll have to post way more code than that to see where the
> >>> mistake is.
> >>>
> >>> John.
> >>>
> >>>> Thanks,
> >>>>
> >>>>
> >>>>
> >>>> _______________________________________________
> >>>> Boost-users mailing list
> >>>> Boost-users_at_[hidden]
> >>>> http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >>> _______________________________________________
> >>> Boost-users mailing list
> >>> Boost-users_at_[hidden]
> >>> http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >>> _______________________________________________
> >>> Boost-users mailing list
> >>> Boost-users_at_[hidden]
> >>> http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >>>
> >> _______________________________________________
> >> Boost-users mailing list
> >> Boost-users_at_[hidden]
> >> http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >> _______________________________________________
> >> Boost-users mailing list
> >> Boost-users_at_[hidden]
> >> http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >>
> > _______________________________________________
> > Boost-users mailing list
> > Boost-users_at_[hidden]
> > http://lists.boost.org/mailman/listinfo.cgi/boost-users
> > _______________________________________________
> > Boost-users mailing list
> > Boost-users_at_[hidden]
> > http://lists.boost.org/mailman/listinfo.cgi/boost-users
> >
> 
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users

> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users