5

Does there exist an approach for calculating the determinant of matrices with low dimensions (about 4), that works well with SIMD (neon, SSE, SSE2)? I am using a hand-expansion formula, which does not work so well. I am using SSE all the way to SSE3 and neon, both under linux. The matrix elements are all floats.

Peter Cordes
  • 328,167
  • 45
  • 605
  • 847
user1095108
  • 14,119
  • 9
  • 58
  • 116
  • 2
    Does this have to be done sequentially or can you calculate multiple independent determinates (e.g. four 4x4 determinate). – Z boson May 04 '15 at 08:53
  • 1
    You also need to be more specific on what determines you want. Do you want exactly 4x4 or do you want 2x2, 3x3, and 4x4? What do you mean by "(about 4)"? – Z boson May 04 '15 at 11:41
  • @Zboson I was looking for a quick way to calculate the determinant, to use with Cramer's rule. I don't have any need for multiple determinants at a time. – user1095108 May 04 '15 at 21:02
  • 1
    Could you explain more of what you're trying to do. What linear systems are you solving? Maybe you can post a little code? – Z boson May 05 '15 at 07:18
  • @Zboson If you check this [page](http://www.euclideanspace.com/maths/algebra/matrix/functions/determinant/fourD/index.htm) you'll see a 4x4 determinant expansion, that baloons into a lot of code. I was simply unnerved that in the time of almost universal SIMD support, I'd use something like that for anything. Looking into my code more closely, it was actually for determining a 4x4 matrix inverse, for some 3D work (I needed the inverse to texture map in software) - there is an algorithm I've seen for the inverse, that does not require calculating a 4x4 determinant at all. – user1095108 May 05 '15 at 08:03
  • 1
    Okay, so we are making progress. So you need a 4x4 matrix inverse and you want to do that with SIMD. But I need a bit more information. Why do you need the inverse? Is it to solve a linear system? Because in that case you may not necessarily need that inverse either. – Z boson May 05 '15 at 08:05
  • @Zboson No it is to do texture mapping of a rectangle. The rectangle is in world coordinates and subject to perspective projection. I look for the inverse of the the projection*view transform to obtain coordinates of textels to map onto the rectangle. Originally, I used Cramer's rule to obtain the inverse, hence the determinant. – user1095108 May 05 '15 at 09:36
  • 2
    In that case I think you should either edit your question (and title) saying you want to determine the inverse of a 4x4 matrix with SSE or make a new question. Also you need to state if this is float or double and also what compiler and OS you are using. – Z boson May 05 '15 at 09:43
  • 2
    Also any other properties of the matrix would be interesting to know. E.g. is it symmetric? Positive definite? Is this an affine transform? – Z boson May 05 '15 at 09:52
  • I mean does your matrix look something like this? https://math.stackexchange.com/questions/65573/inverse-of-a-4x4-matrix#65579 – Z boson May 05 '15 at 09:53
  • @Zboson No, the perspective projective transformation is non-affine, not positive definite either, I think. As for editing the question. I think it's fine, even though I don't need the answers in my code anymore. Basically, I was just looking for a SIMD ways to calculate low-order determinants and I think I've found the answer. If there are others, please contribute (tricks to get the low order determinants efficiently via SIMD). As the answers may span 2 instruction sets, I wasn't looking for specific code. – user1095108 May 05 '15 at 10:30
  • 1
    okay, i understand, if you found an answer to calculate low-order determinants with SIMD please post it. – Z boson May 05 '15 at 11:39

1 Answers1

3

Here's my 5 cents.

determinant of a 2x2 matrix:

that's an exercise for the reader, should be simple to implement

determinant of a 3x3 matrix:

use the scalar triple product. This will require smart cross() and dot() implementations. The recipes for these are widely available.

determinant of a 4x4 matrix:

Use one of the tricks in here. My code:

template <class T>
inline T det(matrix<T, 4, 4> const& m) noexcept
{
  auto const A(make_matrix<T, 2, 2>(m(0, 0), m(0, 1), m(1, 0), m(1, 1)));
  auto const B(make_matrix<T, 2, 2>(m(0, 2), m(0, 3), m(1, 2), m(1, 3)));
  auto const C(make_matrix<T, 2, 2>(m(2, 0), m(2, 1), m(3, 0), m(3, 1)));
  auto const D(make_matrix<T, 2, 2>(m(2, 2), m(2, 3), m(3, 2), m(3, 3)));

  return det(A - B * inv(D) * C) * det(D);
}

determinant of a 5x5+ matrix:

probably use the tricks above.

user1095108
  • 14,119
  • 9
  • 58
  • 116