Inverse of 3x3 Matrix with Complex Numbers
07-15-2019, 11:59 AM
Post: #1 deetee Member Posts: 61 Joined: May 2016
Inverse of 3x3 Matrix with Complex Numbers
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 and stacki). Basic stack operations are done with _add(), _sub() and _mult().
For matrix operations I use matrix A (ma and mai) and B (mb and mbi) and a temporary matrix C (mc and mci).

So far I implemented the following code which works fine:

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

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
07-16-2019, 05:47 AM (This post was last modified: 07-16-2019 05:48 AM by deetee.)
Post: #2 deetee Member Posts: 61 Joined: May 2016
RE: Inverse of 3x3 Matrix with Complex Numbers
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 = ma[r1][c1]; stacki = mai[r1][c1];
stack = ma[r2][c2]; stacki = mai[r2][c2];
stack = ma[r3][c3]; stacki = mai[r3][c3];
_mult(); _mult();
mdet += mult * stack; mdeti += mult * stacki;
}
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 = mdet; stacki = 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 = ma[r1][c1]; stacki = mai[r1][c1]; // Calculates inverse matrix element
stack = ma[r2][c2]; stacki = mai[r2][c2];
stack = ma[r3][c3]; stacki = mai[r3][c3];
stack = ma[r4][c4]; stacki = mai[r4][c4];
_mult(); _rot(); _mult(); _rotup(); _sub(); _push();
stack = mdet; stacki = mdeti;
_div();
mc[r][c] = stack; mci[r][c] = stacki;
}
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
07-16-2019, 10:34 AM
Post: #3
 toml_12953 Senior Member Posts: 1,136 Joined: Dec 2013
RE: Inverse of 3x3 Matrix with Complex Numbers
(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 = ma[r1][c1]; stacki = mai[r1][c1];
stack = ma[r2][c2]; stacki = mai[r2][c2];
stack = ma[r3][c3]; stacki = mai[r3][c3];
_mult(); _mult();
mdet += mult * stack; mdeti += mult * stacki;
}
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 = mdet; stacki = 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 = ma[r1][c1]; stacki = mai[r1][c1]; // Calculates inverse matrix element
stack = ma[r2][c2]; stacki = mai[r2][c2];
stack = ma[r3][c3]; stacki = mai[r3][c3];
stack = ma[r4][c4]; stacki = mai[r4][c4];
_mult(); _rot(); _mult(); _rotup(); _sub(); _push();
stack = mdet; stacki = mdeti;
_div();
mc[r][c] = stack; mci[r][c] = stacki;
}
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.

Tom L

People may say I'm inept but I consider myself to be totally ept.
07-16-2019, 12:07 PM (This post was last modified: 07-16-2019 12:11 PM by deetee.)
Post: #4 deetee Member Posts: 61 Joined: May 2016
RE: Inverse of 3x3 Matrix with Complex Numbers
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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)