Tue 22 Nov 2016 09:57:28 PM UTC, comment #9:
This issue seems to be fixed in Octave 4.2.0 (unless you use a lapack version <3.5.0).
Running the sample code from comment #0 in Octave (on my not so new Ubuntu 14.04 with lapack 3.5.0) gives me an error of only 3.9540e-15. As said in comment #7 this is probably due to my (high enough) lapack version. Older lapack versions probably still have this issue. But 3.5.0 is now (as time has passed by) not very new any more, current version is 3.6.1 and mxe-octave is also built with 3.5.0 today.
I have not checked comment #8.
Is it too early to close this issue? This was always a lapack bug, no Octave bug, when I understand the discussion right. And this bug has been fixed in lapack more than 2.5 years ago, now.
|
Tue 15 Dec 2015 06:07:15 PM UTC, comment #8:
Here is the case where expm from GNU Octave fails to
give the correct results:
A = [0 1e-8 0;-(2e10+4e8/6) -3 2e10;200/3 0 -200/3];
It becomes clear that when we vary the order of magnitude
of the input matrix elements the accuracy of
the results vary.
|
Thu 05 Jun 2014 11:59:24 AM UTC, comment #7:
OK,
I tested LAPACK from 3.4.1 to 3.5.0 (fortran version) with tmpa matrix and the problem is fixed in 3.5.0. So, I presume you have 3.5.0.
About Octave:
*) balance in LAPACK < 3.5.0 has a problem, therefore any Octave version linked to those libraries has the same problem
*) expm works even without balancing. Shall we remove the balancing/unbalancing steps? Shall we make a comment in expm.m?
Marco
|
Thu 05 Jun 2014 10:31:28 AM UTC, comment #6:
I do not see much in octave_config_info (except from -llapack flag that does not help I think).
Running your code gives the following answer:
ans = 0.25000
ans = 0.12500
Thomas
|
Thu 05 Jun 2014 10:00:03 AM UTC, comment #5:
I don't know how to check lapack version (maybe octave_config_info?).
Anyway, can you post the result of
I get
Marco
|
Thu 05 Jun 2014 08:40:02 AM UTC, comment #4:
Hi,
I am not sure how to tell which LAPACK I am using. But I run the test with the MXE for windows (3.8.1-2 installer from http://mxeoctave.osuv.de/) and the norm has the value
3.9589e-015
The norm had more or less the same value when I tested on my debian machine. There, I use the current "testing" package 3.8.1 which I think is using ATLAS.
If you tell me how I can provide more info, I would be happy to.
Best
|
Thu 05 Jun 2014 07:47:46 AM UTC, comment #3:
Hi,
the bug is not solved in 3.8.0. I still see the difference 0.0953. Can you please tell which LAPACK are you using?
Thanks,
Marco
|
Thu 05 Jun 2014 07:34:43 AM UTC, comment #2:
This bug is solved in version 3.8. Could you please close it.
Best
|
Thu 29 Aug 2013 12:08:18 PM UTC, comment #1:
The problem is in the balancing of the matrix.
For instance, the first row of the "balanced" matrix has 1-norm 5e-2 and the first column 9e14 (for the original matrix they are 3.5e-1 and 1e-1).
I think balance uses xGEBAL and that the problem is there (I'm using LAPACK 3.4.1 ang getting similar results with Matlab, where it is written "Balancing can destroy the properties of certain matrices; use it with some care. If a matrix contains small elements that are due to roundoff error [our case], balancing might scale them up to make them as significant as the other elements of the original matrix.").
I see here http://www.netlib.org/lapack/Errata/index2.html that xGEBAL still has some problems in LAPACK 3.4.2.
If you remove the balancing stage, expm is fine. Possible solutions:
1) remove balancing from expm (anyone has a counterexample?)
2) try to fix balance: if, after balancing, the matrix is more unbalanced than the original (for whatever reason), give back the original matrix and a warning.
Marco
|
Tue 27 Aug 2013 01:57:48 PM UTC, original submission:
Hi,
I have ecountered some problem using expm function. The results I obtain have low accuracy.
To reproduce the the problem, I have attached two files:
-expm_data is a data file containing the matrices
tmpa : test matrix
expm_tmpa: result of expm(tmpa) in my PC
expa : expected result (checked both with matlab and an external fortran library)
-expm_test.m is a small script file loading the data and calculating the norm between the expected and the actual result.
In my case the norm has value 0.0953.
best regards
|