Hi all,

Determinant of a matrix can be very very tiny so the epsilon to check it should be really tiny as well.

I use epsilon 1.0e-04f but it's not enough to check determinant of a matrix.

Is it safe to use FLT_EPSILON on this case or it's better to do : if(det == 0.0f) ?

Thanks

What you using the determinant for? How do you know if some method is enough for you?

I got prob on the polar decomp, using 1.0e-04f got problem, so I changed to "== 0.0f" and all worked fine but a floating point needs to be tested using epsilon to be safe.

Maybe on the determinant it's better to not use epsilon or use FLT_EPSILON on this case ?

What algorithm are you using to compute the polar decomposition? There is no generally acceptable epsilon to test two floating points and the epsilon you should use in your case depends on what you are trying to achieve.

This one :

http://research.cs.wisc.edu/graphics/Courses/838-s2002/Papers/polar-decomp.pdf

http://users.soe.ucsc.edu/~pang/160/f98/Gems/GemsIV/polar_decomp/Decompose.c

I don't see any determinant computation in the paper you have linked.. Are you using it to test if the algorithm has converged? Wouldn't it be easier/better to simply compute some matrix norm (like the max norm or the Frobenius norm or some p-norm)?

Why are you writing the algorithm yourself instead of using some linear algebra library? For example, in Lapack or Eigen you may use easily compute the SVD decomposition and then get the polar decomposition from the resulting matrices.

EDIT: How big are your matrices?

if (det==0.0) {do_rank2(Mk, MadjTk, Mk); break;}

Find on the .c you will see this line, it was the line who caused the problem on the animation system.

Why are you writing the algorithm yourself instead of using some linear algebra library?

Polar Decomp of Ken Shoemake works good, it's used on a lot of system and no dependency, the question is not the decomp but the epsilon on a matrix.

#define MATRIX_INVERSE_EPSILON 1e-14
#define MATRIX_EPSILON 1e-6

It really depends on your algorithm and your numeric types. There's lots of algorithms that evaluate determinants. Some algorithms can handle high tolerances for epsilon, like 1e-12. Using floats instead of doubles limits the accuracy of the arithmetic, so high tolerances are probably not going to mean much. That all has to be taken into account. Like you probably know, you never want to just test for exact equivalence to zero unless you're testing for a very specific case.

det = vdot(Mk[0], MadjTk[0]);
if (det==0.0) {do_rank2(Mk, MadjTk, Mk); break;}

In Ken Shoemake's code above, he's just avoiding dividing by zero exception with the "if (det == 0)" check. It seems the algorithm can work for incredibly small numbers, so you don't need anything else other than just avoiding division by zero. Besides, he's using the vdot() function, which as far as I can tell isn't a determinant algorithm, but rather a dot product operation.

Thanks for all answer.

Good catch that it's only to avoid division by zero.