Inverse of 3x3 Matrix with Complex Numbers

[deetee]

Hello!

My next calculator should be able to work with 3x3 matrices with complex numbers (determinant, invert, add, substract, multiply, ...). That for I want to impement fast explicit formulae (found on Wikipedia) rather than some algorithms.

The calculator itself works with two double stacks for real and imaginary parts of numbers (stack[4] and stacki[4]). Basic stack operations are done with _add(), _sub() and _mult().
For matrix operations I use matrix A (ma[3][3] and mai[3][3]) and B (mb[3][3] and mbi[3][3]) and a temporary matrix C (mc[3][3] and mci[3][3]).

So far I implemented the following code which works fine:

Code:

static void _mdet(void) { // MATRIX DET A
  mdet = mdeti = 0.0;
  stack[2] = ma[0][0]; stacki[2] = mai[0][0];
  stack[1] = ma[1][1]; stacki[1] = mai[1][1];
  stack[0] = ma[2][2]; stacki[0] = mai[2][2];
  _mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
  stack[2] = ma[0][1]; stacki[2] = mai[0][1];
  stack[1] = ma[1][2]; stacki[1] = mai[1][2];
  stack[0] = ma[2][0]; stacki[0] = mai[2][0];
  _mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
  stack[2] = ma[0][2]; stacki[2] = mai[0][2];
  stack[1] = ma[1][0]; stacki[1] = mai[1][0];
  stack[0] = ma[2][1]; stacki[0] = mai[2][1];
  _mult(); _mult(); mdet += stack[0]; mdeti += stacki[0];
  stack[2] = ma[2][0]; stacki[2] = mai[2][0];
  stack[1] = ma[1][1]; stacki[1] = mai[1][1];
  stack[0] = ma[0][2]; stacki[0] = mai[0][2];
  _mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
  stack[2] = ma[2][1]; stacki[2] = mai[2][1];
  stack[1] = ma[1][2]; stacki[1] = mai[1][2];
  stack[0] = ma[0][0]; stacki[0] = mai[0][0];
  _mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
  stack[2] = ma[2][2]; stacki[2] = mai[2][2];
  stack[1] = ma[1][0]; stacki[1] = mai[1][0];
  stack[0] = ma[0][1]; stacki[0] = mai[0][1];
  _mult(); _mult(); mdet -= stack[0]; mdeti -= stacki[0];
  stack[0] = mdet; stacki[0] = mdeti;
}
static void _minv(void) { // MATRIX INVERSE A
  _mdet();
  stack[3] = ma[1][1]; stacki[3] = mai[1][1]; // C00
  stack[2] = ma[2][2]; stacki[2] = mai[2][2];
  stack[1] = ma[2][1]; stacki[2] = mai[2][1];
  stack[0] = ma[1][2]; stacki[3] = mai[1][2];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[0][0] = stack[0]; mci[0][0] = stacki[0];
  stack[3] = ma[0][2]; stacki[3] = mai[0][2]; // C01
  stack[2] = ma[2][1]; stacki[2] = mai[2][1];
  stack[1] = ma[0][1]; stacki[2] = mai[0][1];
  stack[0] = ma[2][2]; stacki[3] = mai[2][2];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[0][1] = stack[0]; mci[0][1] = stacki[0];
  stack[3] = ma[0][1]; stacki[3] = mai[0][1]; // C02
  stack[2] = ma[1][2]; stacki[2] = mai[1][2];
  stack[1] = ma[0][2]; stacki[2] = mai[0][1];
  stack[0] = ma[1][1]; stacki[3] = mai[1][2];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[0][2] = stack[0]; mci[0][2] = stacki[0];
  stack[3] = ma[1][2]; stacki[3] = mai[1][2]; // C10
  stack[2] = ma[2][0]; stacki[2] = mai[2][0];
  stack[1] = ma[1][0]; stacki[2] = mai[1][2];
  stack[0] = ma[2][2]; stacki[3] = mai[2][0];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[1][0] = stack[0]; mci[1][0] = stacki[0];
  stack[3] = ma[0][0]; stacki[3] = mai[0][0]; // C11
  stack[2] = ma[2][2]; stacki[2] = mai[2][2];
  stack[1] = ma[0][2]; stacki[2] = mai[0][0];
  stack[0] = ma[2][0]; stacki[3] = mai[2][2];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[1][1] = stack[0]; mci[1][1] = stacki[0];
  stack[3] = ma[0][2]; stacki[3] = mai[0][2]; // C12
  stack[2] = ma[1][0]; stacki[2] = mai[1][0];
  stack[1] = ma[0][0]; stacki[2] = mai[0][0];
  stack[0] = ma[1][2]; stacki[3] = mai[1][2];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[1][2] = stack[0]; mci[1][2] = stacki[0];
  stack[3] = ma[1][0]; stacki[3] = mai[1][0]; // C20
  stack[2] = ma[2][1]; stacki[2] = mai[2][1];
  stack[1] = ma[1][1]; stacki[2] = mai[1][1];
  stack[0] = ma[2][0]; stacki[3] = mai[2][0];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[2][0] = stack[0]; mci[2][0] = stacki[0];
  stack[3] = ma[0][1]; stacki[3] = mai[0][1]; // C21
  stack[2] = ma[2][0]; stacki[2] = mai[2][0];
  stack[1] = ma[0][0]; stacki[2] = mai[0][0];
  stack[0] = ma[2][1]; stacki[3] = mai[2][1];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[2][1] = stack[0]; mci[2][1] = stacki[0];
  stack[3] = ma[0][0]; stacki[3] = mai[0][0]; // C22
  stack[2] = ma[1][1]; stacki[2] = mai[1][1];
  stack[1] = ma[0][1]; stacki[2] = mai[0][1];
  stack[0] = ma[1][0]; stacki[3] = mai[1][0];
  _mult(); _rot(); _mult(); _rotup(); _sub();
  mc[2][2] = stack[0]; mci[2][2] = stacki[0];
  for (byte j = 0; j < 3; j++) for (byte i = 0; i < 3; i++) { // Aij = 1/det * Cij
      stack[0] = mdet; stacki[0] = mdeti;
      stack[1] = mc[j][i]; stacki[1] = mci[j][i];
      _div();
      ma[j][i] = stack[0]; mai[j][i] = stacki[0];
    }
}

Unfortunately this code "eats up" (too) many valuable bytes. I'm sure it can be optimized and shrinked - at least with some loops. But I must resign to manage such loop indices.

Does anyone of you (very smart guys) see how to automate this code in an efficient way?
Any idea and hint is highly welcome.
Regards deetee

[deetee]

Hi!

At least I found some improvement by myself to reduce the code size by summarizing similar code in subroutines - and call them with (a lot of) index parameters:

Code:

static void _mdetrc(byte r1, byte r2, byte r3, byte c1, byte c2, byte c3, int8_t mult) {
  stack[2] = ma[r1][c1]; stacki[2] = mai[r1][c1];
  stack[1] = ma[r2][c2]; stacki[1] = mai[r2][c2];
  stack[0] = ma[r3][c3]; stacki[0] = mai[r3][c3];
  _mult(); _mult();
  mdet += mult * stack[0]; mdeti += mult * stacki[0];
}
static void _mdet(void) { // MATRIX DET A
  mdet = mdeti = 0.0;
  _mdetrc(0, 1, 2, 0, 1, 2, 1);
  _mdetrc(0, 1, 2, 1, 2, 0, 1);
  _mdetrc(0, 1, 2, 2, 0, 1, 1);
  _mdetrc(2, 1, 0, 0, 1, 2, -1);
  _mdetrc(2, 1, 0, 1, 2, 0, -1);
  _mdetrc(2, 1, 0, 2, 0, 1, -1);
  stack[0] = mdet; stacki[0] = mdeti;
}
static void _minvrc(byte r1, byte r2, byte r3, byte r4, byte c1, byte c2, byte c3, byte c4, byte r, byte c) {
  stack[3] = ma[r1][c1]; stacki[3] = mai[r1][c1]; // Calculates inverse matrix element
  stack[2] = ma[r2][c2]; stacki[2] = mai[r2][c2];
  stack[1] = ma[r3][c3]; stacki[1] = mai[r3][c3];
  stack[0] = ma[r4][c4]; stacki[0] = mai[r4][c4];
  _mult(); _rot(); _mult(); _rotup(); _sub(); _push();
  stack[0] = mdet; stacki[0] = mdeti;
  _div();
  mc[r][c] = stack[0]; mci[r][c] = stacki[0];
}
static void _minv(void) { // MATRIX INVERSE A
  _mdet();
  _minvrc(1, 2, 2, 1, 1, 2, 1, 2, 0, 0); // C00
  _minvrc(0, 2, 0, 2, 2, 1, 1, 2, 0, 1); // C01
  _minvrc(0, 1, 0, 1, 1, 2, 2, 1, 0, 2); // C02
  _minvrc(1, 2, 1, 2, 2, 0, 0, 2, 1, 0); // C10
  _minvrc(0, 2, 0, 2, 0, 2, 2, 0, 1, 1); // C11
  _minvrc(0, 1, 0, 1, 2, 0, 0, 2, 1, 2); // C12
  _minvrc(1, 2, 1, 2, 0, 1, 1, 0, 2, 0); // C20
  _minvrc(0, 2, 0, 2, 1, 0, 0, 1, 2, 1); // C21
  _minvrc(0, 1, 0, 1, 0, 1, 1, 0, 2, 2); // C22
  memcpy(ma, mc, 3 * 3 * sizeof(double)); // A = C
  memcpy(mai, mci, 3 * 3 * sizeof(double));
}

Thanks in advance for any hint.
Regards deetee

[toml_12953]

(07-16-2019 05:47 AM)deetee Wrote:  Hi!

At least I found some improvement by myself to reduce the code size by summarizing similar code in subroutines - and call them with (a lot of) index parameters:

Code:

static void _mdetrc(byte r1, byte r2, byte r3, byte c1, byte c2, byte c3, int8_t mult) {
  stack[2] = ma[r1][c1]; stacki[2] = mai[r1][c1];
  stack[1] = ma[r2][c2]; stacki[1] = mai[r2][c2];
  stack[0] = ma[r3][c3]; stacki[0] = mai[r3][c3];
  _mult(); _mult();
  mdet += mult * stack[0]; mdeti += mult * stacki[0];
}
static void _mdet(void) { // MATRIX DET A
  mdet = mdeti = 0.0;
  _mdetrc(0, 1, 2, 0, 1, 2, 1);
  _mdetrc(0, 1, 2, 1, 2, 0, 1);
  _mdetrc(0, 1, 2, 2, 0, 1, 1);
  _mdetrc(2, 1, 0, 0, 1, 2, -1);
  _mdetrc(2, 1, 0, 1, 2, 0, -1);
  _mdetrc(2, 1, 0, 2, 0, 1, -1);
  stack[0] = mdet; stacki[0] = mdeti;
}
static void _minvrc(byte r1, byte r2, byte r3, byte r4, byte c1, byte c2, byte c3, byte c4, byte r, byte c) {
  stack[3] = ma[r1][c1]; stacki[3] = mai[r1][c1]; // Calculates inverse matrix element
  stack[2] = ma[r2][c2]; stacki[2] = mai[r2][c2];
  stack[1] = ma[r3][c3]; stacki[1] = mai[r3][c3];
  stack[0] = ma[r4][c4]; stacki[0] = mai[r4][c4];
  _mult(); _rot(); _mult(); _rotup(); _sub(); _push();
  stack[0] = mdet; stacki[0] = mdeti;
  _div();
  mc[r][c] = stack[0]; mci[r][c] = stacki[0];
}
static void _minv(void) { // MATRIX INVERSE A
  _mdet();
  _minvrc(1, 2, 2, 1, 1, 2, 1, 2, 0, 0); // C00
  _minvrc(0, 2, 0, 2, 2, 1, 1, 2, 0, 1); // C01
  _minvrc(0, 1, 0, 1, 1, 2, 2, 1, 0, 2); // C02
  _minvrc(1, 2, 1, 2, 2, 0, 0, 2, 1, 0); // C10
  _minvrc(0, 2, 0, 2, 0, 2, 2, 0, 1, 1); // C11
  _minvrc(0, 1, 0, 1, 2, 0, 0, 2, 1, 2); // C12
  _minvrc(1, 2, 1, 2, 0, 1, 1, 0, 2, 0); // C20
  _minvrc(0, 2, 0, 2, 1, 0, 0, 1, 2, 1); // C21
  _minvrc(0, 1, 0, 1, 0, 1, 1, 0, 2, 2); // C22
  memcpy(ma, mc, 3 * 3 * sizeof(double)); // A = C
  memcpy(mai, mci, 3 * 3 * sizeof(double));
}

Thanks in advance for any hint.
Regards deetee

Have you benchmarked your routine vs. the built-in routines? I'd be curious to see if there's a significant difference one way or the other.

[deetee]

Hi Tom!

No I didn't benchmark the explicit method and an algorithm based method, but I think the explicit method is much faster. On my Arduino (ATMEGA328P) the result was shown immediately.

On the other hand algorithms are more flexible regarding the (matrix) dimensions.

My intention was to solve matrices with complex numbers up to 3x3. The input of 9 complex numbers is quite painful. For bigger dimensions one should rather use a PC.

Regards
deetee

PS: BTW I'm not using an existing calculator. I try to build my own hardware.