NXC: Matrix algebra library

[HaWe]

[OT]
what is this BranchTest and this BranchComp for? Just for testing? Can it be dropped just to run the code?
(Already ArrayIndex and ArrayReplace I don't quite understand.)
I'll gladly test where I can, but I need to know what exactly I shall be testing. I must admit that I slowly don't quite understand anymore what's happening to my original code.

I want to add 3 aspects (not offending, and not nitpicking, for heaven's sake!) which is important for an "average" NXC user and even for those who understand and try a little more (like maybe me myself) without being an expert or a professional:

1st, any code should still remain "readable" C code,
2nd, of course it's marvellous to have a code that works twice as fast as before; but
3rd, it must not end up in writing pure NBC code to get quite optimized byte code.

To be honest, the fact that nested loops (for (i=...) for (j=....) ...) and array operations (A*=x) and array index operations (A[j]=B[n][m)] are widly unefficient and should be completely replaced by NBC code is not very reassuring.

From the point of view of my current project (stochastic filters) I don't know how important it is to have "maximum speed code" for matrices although it's good to know to have them: I will maybe have to do 3 or 5 or maybe even 10 matrix operations (and for a future particle filter maybe 20 random number generations) every 20msec - the rest are Wait states.
But for a C library it's better to have everything in pure C code (with a little speed disadvantage) than asm code (with a little speed advantage). Too much asm code and tons of #defines will deter NXC users if they don't understand it (e.g., me, in case it's not my own library).

High level code optimization 3+ is good if one knows what it's doing, but e.g., for / while loops or array look ups should be optimized anyway (already by level 0). I personally am suspicious of a code optimization if I don't know what's exactly happening, that's why I'm using code optimization level 1 (only unneeded functions and variables should be erased). The result of any higher optimization level should be accessible by writing pure C code (not asm) without any compiler optimazition, too.

Of course I'm glad that the matrix library leads to this discussion and to the related enhancements of both code and compiler.
Thanks to all for contributing!
[/OT]

[muntoo]

what is this BranchTest and this BranchComp for? Just for testing? Can it be dropped just to run the code?

No!

BranchTest() is similar to:

Code: Select all

if(r != 0)
    goto wherever;

It helps us create an optimized for loop.

And BranchComp() is similar to BranchTest(), but more complicated.

[afanofosc]

These 4 macros that Doc refers to are designed to simply hide the fact that you are essentially writing NBC code and to make the result look exactly like NXC code. No "asm" blocks visible. Just what looks like function calls and labels, which are perfectly valid C code.

BranchTest and BranchComp are both conditional gotos, essentially. You specify the comparison operator (LT, GT, LTEQ, GTEQ, EQ, or NEQ), the label where you should goto, and either one variable (test) or two (compare). Test is "compare to zero" while compare is just a comparison between the two variables. If the comparison or test is true then you goto (aka jump) to the specified label, otherwise you go on to the next statement.

There is nothing at all wrong with the code generated by NXC. The code generated by the very best C compiler will always be slower than hand-written assembly language. I am so very sorry that the code is not perfectly optimized at level 0 or level 1 (i.e., levels that do no optimization at all) and that it is so easy to integrate NBC code into NXC source and to do it in a way that makes it look like plain old ordinary C code. I'm also sorry that people like Doc are afraid to try higher optimization levels and that I am adding new optimizations at higher levels at first so that they can be tested thoroughly before I move them to lower optimization levels.

John Hansen

[muntoo]

I am so very sorry that the code is not perfectly optimized at level 0 or level 1 (i.e., levels that do no optimization at all) and that it is so easy to integrate NBC code into NXC source and to do it in a way that makes it look like plain old ordinary C code.

Is that sarcasm I detect? I'd be happy to try out level 3 optimization on some of my existing programs. (I just hope I'll remember I set it to that!)

[HaWe]

The code generated by the very best C compiler will always be slower than hand-written assembly language.

yes, maybe, although C is reputedly the fastest high level programming language generating insanely fast (machine) code close to asm-generated (machine) code.
Of course I know that NXC is not C, and of course it's fine to be able to integrate asm code to make things a lot faster, but asm is a programming language that many (or even most) NXC users will never use.
As I don't quite understand the latest changes to optimization levels and I don't understand asm code at all (although it's pretty quick):
what code should I take for testing, e.g. using this matrix lib?
The latest asm versions (by John and spillerec) or my old code using just simple C loops and array conversion?
Will the changes affect also my code e.g. of the Mersenne Twister (just doing algebraic bit and byte conversion), the astar (mainly doing array reading/writing), and the FFT (also arithmetics and bit reversals)?

Back to topic, there was a mistake in calculating the Determinant, the Adjugate, and the Inverse for 1x1 matrices.
Those functions use matrices of the kind A[], not of the kind A[][], and as we don't have function overloading the 1x1 dim functions have to be used seperated from n x n dim.
As we don't have recursive functions neither, it's better to use single functions for 2x2 and 3x3 too, because n x n Adjugates can be calculated by the Determinants of the (n-1)x(n-1) sub matrices, which are used for the Transposed Minors resp. Cofactors, and the Transposed Cofactors and the (n x n) Determinate can be used for calculating the Inverse A[n][n]-1.

A[n][n]-1 = det(A[n][n])-1 * Adj(A[n][n]) = det(A[n][n])-1 * (Matrix(Cofactor(A,[n][n])))T

by this it's possible to calculate 3x3 inverse from all 2x2 submatrices, all 4x4 inverse from all 3x3 submatrices, and all 5x5 inverse by all 4x4 submatrices, and so on, pseudo-recursively.

[afanofosc]

While a 1x1 matrix can be represented as a single dimension array it is also perfectly fine to represent it as a 2-d array with 1 row and 1 column. So your code is (or was) just fine wrt 1x1 matrices.

You should use the latest version of the functions that I sent you along with the macros that hide ASM blocks so that your functions all appear to be pure C code. I am working on fourier speedups that I will post later today.

I will upload a new test release today.

John Hansen

[HaWe]

yes, a vector V(x,y) MAY be expanded to a matrix
x , 0
y , 0
but for that you will have to make a type conversion from a 1-dim vector to a 2-dim matrix. THEN you may pass the matrix to a function, and what you get back is again a matrix, which has to be reduced to a 1-dim vector.
That's too complicated.

If your latest compiler version is ready to download, I will test my old code (optimized probably at level 10 by your new compiler) against your asm code.
As I assume, there won't be much difference then in execution speed between both versions any longer. ;)

[afanofosc]

yes, a vector V(x,y) MAY be expanded to a matrix
x , 0
y , 0
but for that you will have to make a type conversion from a 1-dim vector to a 2-dim matrix. THEN you may pass the matrix to a function, and what you get back is again a matrix, which has to be reduced to a 1-dim vector.
That's too complicated.

I do not understand what you mean by x, 0 followed by y, 0. This is not a representation of a 1x1 2d array or at least it doesn't look like one to me. Arrays are not matrices or vectors. They are simply arrays of 1 or more dimensions. A 2d array with 1 row is a vector. It is not a complicated conversion between a vector and a matrix. It's just plain a vector. A 2d array with more than one row and more than one column is a matrix. Using 2-d arrays for both vectors and matrices is a simplification - not a complication. Otherwise you need to have your functions take 2 different types and that needs overloaded functions, etc... which is "too complicated". Just use 2-d arrays where the order is the number of rows in the first dimension. If you ever need a 1d array of one of the rows from a 2-d array it is trivial to get one using ArrayIndex or M[n].

John Hansen

[HaWe]

new functions:
MatrixLinear:

Code: Select all

0  4
1  5 =>  0  1  2  3  4  5  6  7
2  6
3  7

SubMatrix:
deletes row r and column s of a square matrix.

Code: Select all

SubMatrix[2][1] of a 4x4 matrix:
   v
0  4   8   12     
1  5   9   13 
2  6  10   14  <=> [2][1]: start counting at 0; row[2], column[1] 
3  7  11   15
   ^
delete column

deleted column
   v
0       8   12     
1       9   13 
                     <=> deleted row 
3      11   15
 deleted!

0   8   12     
1   9   13    <= new 3x3
3   11  15

testcode:

Code: Select all

#define ArrayInit2D(array, tmp, init_val, dimx, dimy) { \
  ArrayInit(tmp, init_val, dimy);  \
  ArrayInit(array, tmp, dimx);     \
  ArrayInit(tmp,0,0);              \
}


void MatrixLinear(float A[][], float &C[], int R, int S) {
  int i=0, s, r;
  int dim=R*S;
  ArrayInit(C, 0, dim);
  for (s=0; s<S; ++s) {
    for (r=0; r<R; ++r) {
      C[i++]=A[r][s];
    }
  }
}


void SubMatrix(float A[][], float &C[][], int RS, int r, int s) {
  int i=0, j=0, x=0, y=0, dim_1;
  float tmp[];
  
  dim_1=RS-1;
  ArrayInit2D(C, tmp, 0, dim_1, dim_1);
  
  for (y=0; y<RS; ++y) {
    for (x=0; x<RS; ++x) {
      if ((x!=r)&&(y!=s))   {
         C[i++][j]=A[x][y];
         if (i>=dim_1){
           i=0; ++j;
         }
      }
    }
  } 
}
  



task main() {
  float A[2][2], B[2][2], C[2][2], x;
  float O[4][4], T[][], L[];

  O[0][0]=0;   O[0][1]=4;   O[0][2]=8;   O[0][3]=12;
  O[1][0]=1;   O[1][1]=5;   O[1][2]=9;   O[1][3]=13;
  O[2][0]=2;   O[2][1]=6;   O[2][2]=10;  O[2][3]=14;
  O[3][0]=3;   O[3][1]=7;   O[3][2]=11;  O[3][3]=15;

  NumOut(0,56,O[0][0]);  NumOut(24,56,O[0][1]);  NumOut(48,56,O[0][2]);  NumOut(72,56,O[0][3]);
  NumOut(0,48,O[1][0]);  NumOut(24,48,O[1][1]);  NumOut(48,48,O[1][2]);  NumOut(72,48,O[1][3]);
  NumOut(0,40,O[2][0]);  NumOut(24,40,O[2][1]);  NumOut(48,40,O[2][2]);  NumOut(72,40,O[2][3]);
  NumOut(0,32,O[3][0]);  NumOut(24,32,O[3][1]);  NumOut(48,32,O[3][2]);  NumOut(72,32,O[3][3]);

  SubMatrix(O,T,4, 2,1);  // expanded to A[i][j]
  TextOut(0,24,"SubMatrix[2][1]");

  NumOut(0,16,T[0][0]);  NumOut(24,16,T[0][1]);  NumOut(48,16,T[0][2]);
  NumOut(0,08,T[1][0]);  NumOut(24,08,T[1][1]);  NumOut(48, 8,T[1][2]);
  NumOut(0, 0,T[2][0]);  NumOut(24, 0,T[2][1]);  NumOut(48, 0,T[2][2]);

  while (!ButtonPressed(BTNCENTER, false));
  while ( ButtonPressed(BTNCENTER, false));
  ClearScreen();


  while(1);

}

[HaWe]

new functions:

- MatrixMinor[r,s] (= Sub-Determinate[r,s], the determinate of the sub-Matrix[r,s])

- MatrixCofactor[r,s] (MatrixMinor[r,s] * (-1)r+s

- Determinate of a 4x4 Matrix: MatrixDet4x4

Code: Select all

float MatrixMinor(float A[][], int RS, int r, int s) {
   float C[][];

   SubMatrix(A, C, RS, r, s);

   if(RS==3) {
     return MatrixDet2x2(C);
   }
   else
   if(RS==4) {
     return MatrixDet3x3(C);
   }
}


float MatrixCofactor(float A[][], int RS, int r, int s) {
   float x;

   return pow(-1, r+s)* MatrixMinor(A, RS, r, s);;
}


float MatrixDet4x4(float A[][])  {
   float det=0, x=0, M[][];
   int i, j, n0, max0, imax0, jmax0;
   
   for (j=0; j<4; ++j) {    // find column with most 0's: jmax0
     n0=0;
     for (i=0; i<4; ++i) {
       if (A[i][j]==0) ++n0;
     }
     if (max0<n0) {max0=n0; jmax0=j;}
   }
   
   for (i=0; i<4; ++i) {    // find row with most 0's: imax0
     n0=0;
     for (j=0; j<4; ++j) {
       if (A[i][j]==0) ++n0;
     }
     if (max0<n0) {max0=n0; imax0=i;}
   }


   if (jmax0>imax0) {
     for (i=0; i<4; ++i) {         // develop from jmax0
       if ( A[i][jmax0]!=0) {      // 0*?=0
          x=A[i][jmax0]*MatrixCofactor(A, 4, i, jmax0);
          det+=x;
       }
     }
     return det;
   }
   else  {
     for (j=0; j<4; ++j) {         // develop from imax0
       if ( A[imax0][j]!=0) {      // 0*?=0
          x=A[imax0][j]*MatrixCofactor(A, 4, imax0, j);
          det+=x;
       }
     }
     return det;
   }
}

TestCode: see TO post (first untweaked version)