NXC: Matrix algebra library

[HaWe]

both original C code version and tweaked asm code version (some changes switched back for code compatibility) now have been synchronized and can be tested by - more or less - the same task main(): see TO post!
Notice that the last functions according to submatrices and higher level determinates (Minor, Cofactor, 4x4) have NOT been tweaked yet.

To-do list:
4x4 inverse
all functions for 2-dim matrices have to be adjusted and re-written vor 1-dim vectors...

John, when will we have function overloading?

[HaWe]

speed test of both versions (C, asm) with this test code:

Code: Select all

// Matrix algebra / arithmetics
// vers. 1.07a


// Matrix A * scalar: <=>    A*=scalar;

// Matrix A + Matrix B: <=>  A+= B;


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


void MatrixMatrixMult(float A[][], float B[][], float &C[][], int N, int M, int K){
  int i, j, s;
  float tmp[];
  ArrayInit2D(C,tmp, 0, N, K);
  for (i=0; i<N; ++i) {
    for (j=0; j<K; ++j) {
       C[i][j]=0;
       for (s=0; s<M; ++s) {
         C[i][j]=C[i][j] + A[i][s]*B[s][j];
      }
    }
  }
}

                                                // Transposed / Transponierte
void MatrixTransp(float A[][], float &C[][], int R, int S) {
   int r, s;
   float tmp[];
   ArrayInit2D(C,tmp, 0, S,R);

   for (s=0; s<S; s++) {
     for (r=0; r<R; r++) {
       C[s][r]=A[r][s];
     }
   }
}




void MatrixRotate2DMath(float A[][], float &C[][], float angleDeg) {
   float tmp[], x, B[][];

   ArrayInit2D(C, tmp, 0, 2, 2);

   x=cosd(angleDeg);
   ArrayInit2D(B, tmp, x, 2, 2);
   x=sind(angleDeg);
   B[0][1]=-x;
   B[1][0]= x;

   MatrixMatrixMult(B, A, C, 2, 2, 2);
}


#define MatrixRotate2DGeo(A,C,alpha)  MatrixRotate2DMath(A,C,-alpha)



float MatrixDet1x1(float A[]) {                 // Determinante 1x1
   return ( A[0]);
}


float MatrixDet2x2(float A[][]) {               // Determinante 2x2
   return ( A[0][0]*A[1][1]- A[0][1]*A[1][0] );
}


float MatrixDet3x3(float A[][]) {               // Determinante 3x3
   return (A[0][0]*A[1][1]*A[2][2]
          +A[0][1]*A[1][2]*A[2][0]
          +A[0][2]*A[1][0]*A[2][1]
          -A[0][2]*A[1][1]*A[2][0]
          -A[0][1]*A[1][0]*A[2][2]
          -A[0][0]*A[1][2]*A[2][1]);
}


void MatrixAdj1x1 (float A[], float &C[]) {     // Adjugate=Adjunkte 1x1
  ArrayInit(C, 1, 1);
}


void MatrixAdj2x2 (float A[][], float &C[][]) {  // Adjugate=Adjunkte 2x2
  float tmp[];
  ArrayInit2D(C, tmp, 0, 2,2);
  C[0][0]=  A[1][1];   C[1][0]= -A[1][0];
  C[0][1]= -A[0][1];   C[1][1]=  A[0][0];
}


void MatrixAdj3x3 (float A[][], float &C[][]) {  // Adjugate=Adjunkte 3x3
   float tmp[];
   ArrayInit2D(C, tmp, 0, 3,3);
   C[0][0]=A[1][1]*A[2][2]-A[1][2]*A[2][1];  C[0][1]=A[0][2]*A[2][1]-A[0][1]*A[2][2];  C[0][2]=A[0][1]*A[1][2]-A[0][2]*A[1][1];

   C[1][0]=A[1][2]*A[2][0]-A[1][0]*A[2][2];  C[1][1]=A[0][0]*A[2][2]-A[0][2]*A[2][0];  C[1][2]=A[0][2]*A[1][0]-A[0][0]*A[1][2];

   C[2][0]=A[1][0]*A[2][1]-A[1][1]*A[2][0];  C[2][1]=A[0][1]*A[2][0]-A[0][0]*A[2][1];  C[2][2]=A[0][0]*A[1][1]-A[0][1]*A[1][0];
}


bool MatrixInverse1x1 (float A[], float &C[]) {
   ArrayInit(C, 0, 1);
   if (C[0]==0) return 0;
   else {
     C[0]=1/A[0];
     return 1;
   }
}

bool MatrixInverse2x2 (float A[][], float &C[][]) {
   float det;
   float tmp[];
   ArrayInit2D(C, tmp, 0, 2,2);

   det=MatrixDet2x2(A);
   if (det==0) return 0;
   else {
     MatrixAdj2x2(A, C);
     C/=det;
     return 1;
   }
}


bool MatrixInverse3x3 (float A[][], float &C[][]) {
   float det;
   float tmp[];
   ArrayInit2D(C, tmp, 0, 3,3);

   det=MatrixDet3x3(A);
   if (det==0) return 0;
   else {
     MatrixAdj3x3(A, C);
     C/=det;
     return 1;
   }
}

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];
         ++i;
         if (i>=dim_1){
           i=0; ++j;
         }
      }
    }
  }
}

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;
   int i, j=0, n0, max0=0, imax0, jmax0=0;

   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;
       }
     }
     //TextOut(72,56,"j"); NumOut(80,56,jmax0); //debug
     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;
       }
     }
     //TextOut(72,56,"i"); NumOut(80,56,imax0); //debug
     return det;
   }
}

// speed test
long time0, runtime;

task main() {

  float A[][], B[][], C[][], x;
  float O[][], T[][], L[][], tmp[];;
  int loop;
  time0=CurrentTick();
  
 for (loop=0; loop<100; ++loop)
 {
  ArrayInit2D(A, tmp, 0, 2,2);
  ArrayInit2D(B, tmp, 0, 2,2);
  ArrayInit2D(C, tmp, 0, 2,2);

  ArrayInit2D(O, tmp, 0, 4,4);

  NumOut(0,56,loop);

  O[0][0]=0;   O[0][1]=0;   O[0][2]=8;   O[0][3]=22;
  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]=22;  O[3][3]=15;

  x=MatrixDet4x4(O);

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


  x=MatrixDet3x3(T);


  SubMatrix(T,C,3, 0,1);

  x=MatrixDet2x2(C);


  x=MatrixMinor(T, 3, 0,1);

  x=MatrixCofactor(T, 3, 0,1);


  A[0][0]=1;   A[0][1]=3;
  A[1][0]=2;   A[1][1]=4;

  x=10;
  B=A; B*=x;

  C=A; C+=B;

  MatrixAdj2x2(A,C);

  MatrixInverse2x2(A,C);

  MatrixInverse2x2(A,B);

  MatrixMatrixMult(A,B,C,2,2,2);


  A[0][0]=1;   A[0][1]=1;
  A[1][0]=1;   A[1][1]=0;


  MatrixRotate2DGeo(A,B,30);

  MatrixRotate2DGeo(A,C,90);

 }
 runtime=CurrentTick()-time0;
 TextOut(0,40,"run time");
 NumOut (0,32, runtime);   TextOut(72,32,"msec");
// 10514 107TEST BCC20110707

 while(1);
}

asm version (updated #defines)

Code: Select all

            // Matrix algebra / arithmetics
// vers. 0.92t (HaWe+jch+splrc version)


#define MatrixMulS(_Out, _A, _s) asm { mul _Out, _A, _s }
#define MatrixDivS(_Out, _A, _s) asm { div _Out, _A, _s }

#define MatrixAdd(_Out, _A, _B) asm { add _Out, _A, _B }
#define MatrixSub(_Out, _A, _B) asm { sub _Out, _A, _B }




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


void MatrixTransp(float A[][], float &C[][], int R, int S) {
   float tmp[], arr_temp[], val_temp;
   int s, r;
   ArrayInit(tmp, 0, R);
   ArrayInit(C, tmp, S);
   s = S;
   lbl_Trans_start_s:
   {
      s--;
      r = R;
      lbl_Trans_start_r:
      {
         r--;
         ArrayIndex(arr_temp, A, r);
         ArrayIndex(val_temp, arr_temp, s);
         ArrayReplace(tmp, r, val_temp);
      }
      BranchTest(GT, lbl_Trans_start_r, r);
      ArrayReplace(C, s, tmp);
   }
   BranchTest(GT, lbl_Trans_start_s, s);
}

                                               // N x M Matrix  *  M x K Matrix, OLD version
void MatrixMatrixMult_OLD(float A[][], float B[][], float &C[][], int N, int M, int K){
  int i, j, s;

  for (i=0; i<N; ++i) {
    for (j=0; j<K; ++j) {
       C[i][j]=0;
       for (s=0; s<M; ++s) {
         C[i][j]=C[i][j] + A[i][s]*B[s][j];
      }
    }
  }
}

                                               // N x M Matrix  *  M x K Matrix, NEW version
void MatrixMatrixMult(float A[][], float B[][], float &C[][], int N, int M, int K){
   float arr_temp[];
   //Set size of C
   ArrayInit( arr_temp, 0, K );
   ArrayInit( C, arr_temp, N );

lbl_Mult_start_i:   N--;   //Loop start
      //Speed up arrays
      float A_i[];
      float C_i[];
      ArrayIndex(A_i, A, N);
      ArrayIndex(C_i, C, N);

      int j = K;
lbl_Mult_start_j: j--;   //Second loop start
         float sum = 0;

         int s = M;
lbl_Mult_start_s: s--;   //Third loop start

            //sum += A_i[s]*B[s][j];
            float temp, sum_temp;
            ArrayIndex(sum_temp, A_i, s);
            ArrayIndex(arr_temp, B, s);
            ArrayIndex(temp, arr_temp, j);
            sum_temp *= temp;
            sum += sum_temp;
         asm { brtst GT, lbl_Mult_start_s, s }
         ArrayReplace(C_i, j, sum);
      asm { brtst GT, lbl_Mult_start_j, j }
      ArrayReplace(C, N, C_i);
     asm { brtst GT, lbl_Mult_start_i, N }
}



void MatrixRotate2DMath(float A[][], float &C[][], float angleDeg) {
   float tmp[], x, B[][];

   ArrayInit2D(C, tmp, 0, 2, 2);

   x=cosd(angleDeg);
   ArrayInit2D(B, tmp, x, 2, 2);
   x=sind(angleDeg);
   B[0][1]=-x;
   B[1][0]= x;

   MatrixMatrixMult(B, A, C, 2, 2, 2);
}


#define MatrixRotate2DGeo(A,C, d) MatrixRotate2DMath(A,C,-d)



float MatrixDet2x2(float A[][]) {               // Determinante 1x1, 2x2, 3x3
   float r0[], r1[], r2[], v00, v01, v02, v10, v11, v12, v20, v21, v22;

      r0 = A[0];
      r1 = A[1];
      return ( r0[0]*r1[1]- r0[1]*r1[0] );
}


float MatrixDet3x3(float A[][]) {               // Determinante  3x3
   float r0[], r1[], r2[], v00, v01, v02, v10, v11, v12, v20, v21, v22;

      r0 = A[0];
      r1 = A[1];
      r2 = A[2];
      v00 = r0[0]; v01 = r0[1]; v02 = r0[2];
      v10 = r1[0]; v11 = r1[1]; v12 = r1[2];
      v20 = r2[0]; v21 = r2[1]; v22 = r2[2];
      return (v00*v11*v22
             +v01*v12*v20
             +v02*v10*v21
             -v02*v11*v20
             -v01*v10*v22
             -v00*v12*v21);

}


void MatrixAdj2x2 (float A[][], float &C[][]) {  // Adjugate=Adjunkte 1x1, 2x2, 3x3
  float r0[], r1[], r2[], c0[], c1[], c2[], t0, t1, t2;

    ArrayIndex(r0, A, 0);
    ArrayIndex(r1, A, 1);
    ArrayIndex(t0, r1, 1);
    ArrayIndex(t1, r0, 1);
    asm { mul t1, t1, -1 }
    ArrayBuild(c0, t0, t1);
    ArrayIndex(t0, r1, 0);
    asm { mul t0, t0, -1 }
    ArrayIndex(t1, r0, 0);
    ArrayBuild(c1, t0, t1);
    ArrayBuild(C, c0, c1);

}


void MatrixAdj3x3 (float A[][], float &C[][]) {  // Adjugate=Adjunkte 1x1, 2x2, 3x3
  float r0[], r1[], r2[], c0[], c1[], c2[], t0, t1, t2;


    float v00, v01, v02, v10, v11, v12, v20, v21, v22;
    r0 = A[0]; r1 = A[1]; r2 = A[2];
    v00 = r0[0]; v01 = r0[1]; v02 = r0[2];
    v10 = r1[0]; v11 = r1[1]; v12 = r1[2];
    v20 = r2[0]; v21 = r2[1]; v22 = r2[2];

    t0 = v11*v22-v12*v21;
    t1 = v02*v21-v01*v22;
    t2 = v01*v12-v02*v11;
    ArrayBuild(c0, t0, t1, t2);

    t0 = v12*v20-v10*v22;
    t1 = v00*v22-v02*v20;
    t2 = v02*v10-v00*v12;
    ArrayBuild(c1, t0, t1, t2);

    t0 = v10*v21-v11*v20;
    t1 = v01*v20-v00*v21;
    t2 = v00*v11-v01*v10;
    ArrayBuild(c2, t0, t1, t2);

    ArrayBuild(C, c0, c1, c2);

}


bool MatrixInverse2x2(float A[][], float &C[][]) {
   float det;

   det=MatrixDet2x2(A);
   if (det==0) return 0;
   else {
     MatrixAdj2x2(A, C);
     C/=det;
     return 1;
   }
}


bool MatrixInverse3x3(float A[][], float &C[][]) {
   float det;

   det=MatrixDet3x3(A);
   if (det==0) return 0;
   else {
     MatrixAdj3x3(A, C);
     C/=det;
     return 1;
   }
}



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];
         ++i;
         if (i>=dim_1){
           i=0; ++j;
         }
      }
    }
  }
}


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;
   int i, j=0, n0, max0=0, imax0, jmax0=0;

   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;
       }
     }
     //TextOut(72,56,"j"); NumOut(80,56,jmax0); //debug
     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;
       }
     }
     //TextOut(72,56,"i"); NumOut(80,56,imax0); //debug
     return det;
   }
}


// speed test
long time0, runtime;

task main() {

  float A[][], B[][], C[][], x;
  float O[][], T[][], L[][], tmp[];;
  int loop;
  time0=CurrentTick();

 for (loop=0; loop<100; ++loop)
 {
  ArrayInit2D(A, tmp, 0, 2,2);
  ArrayInit2D(B, tmp, 0, 2,2);
  ArrayInit2D(C, tmp, 0, 2,2);

  ArrayInit2D(O, tmp, 0, 4,4);

  NumOut(0,56,loop);

  O[0][0]=0;   O[0][1]=0;   O[0][2]=8;   O[0][3]=22;
  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]=22;  O[3][3]=15;

  x=MatrixDet4x4(O);

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


  x=MatrixDet3x3(T);


  SubMatrix(T,C,3, 0,1);

  x=MatrixDet2x2(C);


  x=MatrixMinor(T, 3, 0,1);

  x=MatrixCofactor(T, 3, 0,1);


  A[0][0]=1;   A[0][1]=3;
  A[1][0]=2;   A[1][1]=4;

  x=10;
  B=A; B*=x;

  C=A; C+=B;

  MatrixAdj2x2(A,C);

  MatrixInverse2x2(A,C);

  MatrixInverse2x2(A,B);

  MatrixMatrixMult(A,B,C,2,2,2);


  A[0][0]=1;   A[0][1]=1;
  A[1][0]=1;   A[1][1]=0;


  MatrixRotate2DGeo(A,B,30);

  MatrixRotate2DGeo(A,C,90);

 }
 runtime=CurrentTick()-time0;
 TextOut(0,40,"run time");
 NumOut (0,32, runtime);   TextOut(72,32,"msec");
// 10514 107TEST BCC20110707

 while(1);
}

result (updated)
(c-like 107TEST BCC20110707 )
opt.lev 1: 10514 mesc
opt.lev 2: 9158 mesc
opt lev. 6: 9155 mesc

(asm 092TEST BCC20110707)
opt.lev 1: 8464 msec
opt.lev 2: 7388 msec
opt lev. 6: 7384 mesc

(c-like 107TEST BCC20110719 )
opt.lev 1: 10522 msec
opt.lev 2: 9038
opt lev. 3: 8061
opt.lev 4: 8061
opt.lev 5: 8061
opt lev. 6: 8061


(asm 092TEST BCC20110719)
opt.lev 1: 8464 msec
opt.lev 2: file error -1
opt lev. 3: file error -1
opt.lev 4: file error -1
opt.lev 5: file error -1
opt lev. 6: file error -1

[afanofosc]

Speed tests should not include large amounts of unrelated code and you should never have LCD drawing involved. The unrelated code and especially drawing to the LCD will invalidate your results.

John Hansen

[HaWe]

oh, I though NumOut just writes to a buffer?
but anyhow, over all it writes 100 times to the display, and unrelated code is not much (unused code is supposed to be erased though) - I think that's a quite normal test for normal programs under normal conditions... :)

[HaWe]

ps:
without NumOut it's 7373 ms for the asm version (vs. 7384) - that's really not much enough to talk about...^^

[afanofosc]

all functions for 2-dim matrices have to be adjusted and re-written vor 1-dim vectors...

Functions designed to work on matrices of any size (1x1, 2x2, 3x3, etc...) will always take two dimensional arrays since a matrix has two dimensions (i.e., row = dimension 1 and column = dimension 2). Functions designed to work with vectors are not the same beast as functions designed to work with matrices. A matrix with 1 row and N columns or N rows and 1 column is still a matrix and can easily and simply be represented as a two dimensional array.

Why do you insist that a 1x1 matrix should be represented using a one dimension array? It doesn't make any sense to me.

John Hansen

[HaWe]

John,
what I mean is:
matrix and vector actually are the same.
I'm no mathematics teacher, maybe a mathematician can eplain it better.

I'm sorry that I maybe not can discribe this sufficiently enough with my instruments of English.

of course I have to rewrite those new functions anyway, but then I have to call them differently than the related 2-dim ones, because of no function overloading. That's the only "problem".

[linusa]

I think John is saying that you can always use a 2x1 matrix as column vector, and a 1x2 matrix as row vector.

So when your 2D rotation operation is with
and
,

you always have 2-dimensional arrays -- for x and y, too. Even if 1 dimension only has 1 element. And a scalar is a 2dim-array with just 1 element each. That way you don't need overloading at all.

If performance is a problem in your application, you optimize at the very last moment.

Edit: removed comment about the post being in the wrong thread -- thanks for moving

[HaWe]

what you mean with your set formulas and R2x2 is beyond me.
??
For NXC
A[n][1] is != A[n]
if a function is designed for
void foo(int A[][]), the compiler does accept B[6][4] and C[6][1], but does not accept D[6].

On the other hand, I don't want to expand all 1-dim arrays[n] to arrays[n][1] before I can use them and back to calculate array[n']
(In C it's possible by using int * A) .

Of course it's possible to do that, of course I can use temporary variables, of course I can reproduce all formulas to apply to 1-dim arrays, and of course I can use different other workarounds too, but my wish is that I may use any function directly without intermediate temp variables and without having to write spaghetti code.

Of course I know, too, that my wish currently is unattainable - but it's just a wish.

[linusa]

I noticed that the whole discussion is moving towards the direction "how to implement a linear algebra math library", and how to "properly support NxM matrices with high flexibility and good performance."

Well, OpenCV has done an excellent job there, for example. With both C and C++ datatypes. I suggest looking at the OpenCV book, it explains the C datastructure for matrices and images (in OpenCV, images are just big matrices). It also explains the reasons WHY this was implemented the way it is. The book (title "Learning OpenCV", O'Reilly) only focuses on OpenCV 1.3 and its C-interface.

Anyway, the solution is to handle this problem object-oriented (like it or not, doc-helmut). You design a basic matrix data-type, like

Code: Select all

struct MyFloatMatrix
  float data[][];
  int M;
  int N;
  bool transposed;
};

You always pass this datatype to functions. They immediately know the size (MxN) of the matrix, so they can also do sanity-checks and assertions in debug mode, if you can for example multiply 2 matrices. I built in the transposed-flag I talked about.

If you want an easy way to get sub-matrices, you add more info to this struct:

Code: Select all

struct MyFloatMatrix
  float data[][];
  int M;
  int N;
  bool transposed;
  int ROI_startRow;
  int ROI_startCol;
  int ROI_numRows;
  int ROI_numCols;
};

If somebody wants the top-left 3x2 matrix from a 4x4 one, you simply set the ROI (region of interest, the term used in OpenCV) params. In this case, you'd set ROI_startRow to 0, ROI_startcol to 0, ROI_numRows to 3, and ROI_numCols to 2.

All functions handling matrix data interpret this ROI settings. For example, to access all matrix elements, you'd do like this

Code: Select all

for (int i = M.ROI_startRow; i < M.ROI_numRows; ++i) {
  for (int j = M.ROI_startCol; j < M.ROI_numCols; ++j) {
     // do stuff to i, j, which are the correct indices for the submatrix
  }//end for
}//end for

Another big advantage (not for NXC though) is to use pointers for the data-array and to use lazy copying: If you have submatrix, it still referes to the original matrix-data, with just new ROI settings.

This makes operations that need submatrices very cheap (almost for free, CPU and memory wise).

OpenCV also implements some hardcoded "3x3" etc matrix data types. But then they have LOTs of individual functions. They are probably generated.

This is the trade-off. You either have very flexible code, which might have some overhead, or you do hardcoded MxN matrix data types, and have lots of generated or copy-pasted functions.
Edit: removed comment about the post being in the wrong thread -- thanks for moving