[Insight-developers] SVD compiled with optimziation
Miller, James V (Research)
millerjv at crd.ge.com
Tue Nov 29 10:10:18 EST 2005
Karthik,
What we have done in the past for slamch.c and dlamch.c is place
pragmas in the code to turn off the optimizations. This is done
for the Microsoft and Intel compilers.
Does dsvdc.c end up calling dlamch (or a similar routine to determine
machine precision)?
The linpack/lapack/eispack creators were VERY careful with machine
precision. Much more so than we are in image analysis. However, processor
advances confuse some of this old code. For instance, the Intel processors
have extra wide floating point registers. Wider than IEEE. So whenever
a calculation can be performed and stored in registers the result will be
of higher precision than if any of the temporaries were stored in memory.
This old code usually performs tests to determine the operating precision
of the machine to determine tolerances. There is really no telling whether
the precision determination will be done within registers or within memory.
The optimization level of the compiler will change how registers are used
and how often temporary results are stored in registers instead of in memory.
I can appreciate such behavior causing a routine like svd to completely
fail. Recall that the purpose of svd is to identify the singular
values.
Jim
-----Original Message-----
From: insight-developers-bounces+millerjv=crd.ge.com at itk.org
[mailto:insight-developers-bounces+millerjv=crd.ge.com at itk.org]On Behalf
Of Karthik Krishnan
Sent: Monday, November 28, 2005 8:30 PM
To: Insight Developers List
Cc: suyash at cs.utah.edu
Subject: [Insight-developers] SVD compiled with optimziation
Hello,
I ran into the same problem today, using itk::ImageMomemtsCalculator. (I'm using gcc40 on debian)
Turns out, that uses vnl_svd which uses netlib/dsvdc.c and all these are related.
Funnily enough, turns out that is some numerical issue in netlib/dsvdc when optimization is turned on.
Did you compile with the -O3 flag (usually the Release mode on CMake) ? For now, I changed the Utilities/vxl/netlib/Makefile to -O2 on the build line that compiles dsvdc.c
and everything was working great.
I haven't snooped around to see what exactly is causing the problem, but there seem to be others in the VXL community who've run into similar issues.
The following snippet in dsvdc.c and the cvs logs from the vxl sources.
<snip>
----------------------------
revision 1.6
date: 2002/03/13 15:51:28; author: fsm; state: Exp; lines: +12 -2
hack to avoid excess precision
----------------------------
/*
* Calling this ensures that the operands are spilled to
* memory and thus avoids excessive precision when compiling
* for x86 with heavy optimization (gcc). It is better to do
* this than to turn on -ffloat-store.
*/
static int fsm_ieee_doubles_differ(double *x, double *y)
{
return *x != *y;
}
</snip>
In any case comparing two doubles for equality sounds questionable.
What's bothersome is that the inverse with -O3 builds is not marginally wrong, but quite wrong for a 3x3 rotation matrix ?
Has anybody run into this issue. Is it limited only to gcc ?
I've changed my nightly builds to reflect this. All the associated (23 tests or so fail now on the release builds)
http://www.itk.org/Testing/Sites/Sabbath.kitware/zRel22-Linux-gcc40/20051128-0100-Nightly/Test.html
Thanks
karthik
---------------------------------------------
Suyash,
You might not have enough precision in the matrix values that you are reading from the command line. Try generating the values with more
precision from whatever tool is generating them. Try going out to about 12 decimal places.
I have seen this before where one program computes a matrix that is supposed to be invertible, writes the matrix in ascii, another program reads that matrix and it is no longer invertible.
Jim
-----Original Message-----
From: insight-users-bounces+millerjv=crd.ge.com at itk.org <http://www.itk.org/mailman/listinfo/insight-users> [mailto:insight-users-bounces+millerjv=crd.ge.com at itk.org <http://www.itk.org/mailman/listinfo/insight-users>]On Behalf Of Suyash P. Awate
Sent: Friday, November 11, 2005 1:06 PM
To: jiafucang; jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users>; insight-users at itk.org <http://www.itk.org/mailman/listinfo/insight-users>
Subject: [Insight-users] Re: On incorrect matrix inverse
Hi everybody,
Thanks for the generous help.
I was also trying to get the inverse concerning with the registration application, like Fucang (see email attached).
Here is my piece of code that gave me the wrong result:
------------
const unsigned int Dimension = 3;
typedef double PixelType;
typedef itk::Matrix <PixelType, Dimension, Dimension> MatrixType;
MatrixType matrix;
matrix.Fill (0);
matrix (0,0)= atof (argv[4]);
matrix (0,1)= atof (argv[5]);
matrix (0,2)= atof (argv[6]);
matrix (1,0)= atof (argv[7]);
matrix (1,1)= atof (argv[8]);
matrix (1,2)= atof (argv[9]);
matrix (2,0)= atof (argv[10]);
matrix (2,1)= atof (argv[11]);
matrix (2,2)= atof (argv[12]);
vnl_matrix_fixed <PixelType, Dimension, Dimension> inverse;
inverse= matrix.GetInverse();
-----------
matrix
0.99925 0.0371382 -0.0109982
-0.0373772 0.99905 -0.0223916
0.0101562 0.0227859 0.999689
warnings:
InsightToolkit-2.2.0/Utilities/vxl/core/vnl/algo/vnl_svd.txx: suspicious return value (2) from SVDC
InsightToolkit-2.2.0/Utilities/vxl/core/vnl/algo/vnl_svd.txx: M is 3x3
M = [ ...
0.9992500000000 0.0371382000000 -0.0109982000000
-0.0373772000000 0.9990500000000 -0.0223916000000
0.0101562000000 0.0227859000000 0.9996890000000 ]
matrix inverse using vnl in itk (incorrect)
-0.982557 -0.0644525 -0.174434
-0.0677511 0.997618 0.0130162
-0.17318 -0.0246073 0.984583
matrix inverse in matlab (correct)
0.9992 -0.0374 0.0102
0.0371 0.9991 0.0228
-0.0110 -0.0224 0.9997
------------
I had also tried using the AffineTransform but for some reason could not get that to work and so I wrote this simple matrix code. Perhaps that was due to the same problem that Fucang reported.
Currently, however, my problem is different. As we can see I am not using any Transform class at all. Moreover I am getting warnings from the vnl_svd class, so vnl is definitely involved here.
Julien: I am giving 9 numbers (truncated: 0.99925 0.0371382 -0.0109982 -0.0373772 0.99905 -0.0223916 0.0101562 0.0227859 0.999689) to the matrix in itk and also the same numbers to the matrix in matlab. I got these numbers from the output of a separate program.
Just FYI, I have compiled ITK on Suse 10.0 using GCC 4.0. I wonder if that has anything to do with this !
Best regards,
Suyash.
----- Original Message -----
From: <mailto:jiafucang at asisz.com <http://www.itk.org/mailman/listinfo/insight-users>> jiafucang
To: <mailto:suyash at cs.utah.edu <http://www.itk.org/mailman/listinfo/insight-users>> suyash at cs.utah.edu <http://www.itk.org/mailman/listinfo/insight-users> ; <mailto:jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users>> jjomier at cs.unc.edu <http://www.itk.org/mailman/listinfo/insight-users>
Sent: Friday, November 11, 2005 7:10 AM
Subject: On incorrect matrix inverse
Hello, Suyash, Jomier,
I have reported similar problem on
<http://public.kitware.com/pipermail/insight-users/2005-October/015441.html> http://public.kitware.com/pipermail/insight-users/2005-October/015441.html
and added a bug 2450
<http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24 <http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24>> http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24 <http://www.itk.org/Bug/bug.php?op=show&bugid=2450&pos=24>
I think it is because ITK does not computer inverse when setparameters used, and not VNL problem.
HTH
Fucang
_______________________________________________
Insight-developers mailing list
Insight-developers at itk.org
http://www.itk.org/mailman/listinfo/insight-developers
More information about the Insight-developers
mailing list