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@lists.boost.org> on behalf of John Maddock <jz.maddock@googlemail.com> Sent: Wednesday, May 4, 2016 4:49:07 PM To: boost-users@lists.boost.org 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@lists.boost.org> on behalf of John Maddock <jz.maddock@googlemail.com> Sent: Wednesday, May 4, 2016 12:29:31 PM To: boost-users@lists.boost.org 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@lists.boost.org> on behalf of John Maddock <jz.maddock@googlemail.com> Sent: Tuesday, May 3, 2016 6:27:39 PM To: boost-users@lists.boost.org Subject: Re: [Boost-users] converting to float128 from double
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
On 03/05/2016 17:35, Cooper, Bridgette R D wrote: 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@lists.boost.org> on behalf of John Maddock <jz.maddock@googlemail.com> Sent: Tuesday, May 3, 2016 5:06:47 PM To: boost-users@lists.boost.org 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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
_______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users