[Insight-developers] SVD compiled with optimziation
Karthik Krishnan
Karthik.Krishnan at kitware.com
Tue Nov 29 13:05:52 EST 2005
Miller, James V (Research) wrote:
>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)?
>
>
>
I don't know if they called dlamch etc., but I don't think so.
I ended up manually editing the netlib/Makefile to add
-fno-inline-functions to the already existing flags
-w -O3 -fPIC
and that took care of the issue.
Since I added these flags to just the files dsvdc.o, csvdc.o, zsvdc.o,
ssvdc.o and not to any of the others (dlamch.o or machineparams.o), I am
guessing that may not be the reason. I still think its just that the
function
static int fsm_ieee_doubles_equal(const doublereal *x, const doublereal *y);
[ which as you mentioned seems to be defined to explicitly typecast the
high precision registers to doubles is being inlined away with -O3 and
that defeats the purpose of this workaround function ]. Is it possible
that woraround only with gcc40, which got smart enough to figure out the
workaround.
from gcc docs
|-O3|
Optimize yet more. -O3 turns on all optimizations specified by -O2
and also turns on the -finline-functions, -funswitch-loops and
-fgcse-after-reload options.
Thanks
Karthik
>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
>_______________________________________________
>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