EV3 CSLite C libs: matrix algebra, array sort and more

[gloomyandy]

Doc,
most programmers would not choose to do matrix manipulations in C (they would use Java, python, C++ or some other higher level language), and for those that do, the code I posted that shows how to do it using pointers is simple code...

[HaWe]

Maybe I have overlooked something but I didn't find something as simple as

Code: Select all

void MatrixMatrixMult(float A[][], float B[][], float &C[][], int N, int M, int K)

yet which simply computes the matrix product AxB=C.

I once started programming by Turbo Pascal which generated executables for MSDOS - no VM at all .
Turbo Pascal and Delphi also had an API for Fischertechnik robotic computing, I didn't have to struggle with makefiles and project managers and .so and .o, I just pressed a key and it compiled and it linked unseen and it immediately worked (and also Turbo Pascal was featuring pointers, I wrote my own quick and shell sort programs, e.g. for doubly linked lists ).
I could write programs for robots and didn't have to invent a programming language first.

- IIRC, life was so much easier in those days.

[gloomyandy]

Doc,
the code in the following link:
http://www.mymathlib.com/c_source/matri ... matrices.c
Which I posted above. Does pretty much exactly what your matrix multiply does (the parameters are in a different order but I'm sure you can cope with that). There is an example that shows how it can be used. Did you look at the code? What is your problem with it?

[HaWe]

Andy,
thanks a lot for your help - indeed this is the code I have been looking for.

Code: Select all

////////////////////////////////////////////////////////////////////////////////
// File: multiply_matrices.c                                                  //
// Routine(s):                                                                //
//    Multiply_Matrices                                                       //
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  void Multiply_Matrices(double *C, double *A, int nrows, int ncols,        //
//                                                    double *B, int mcols)   //
//                                                                            //
//  Description:                                                              //
//     Post multiply the nrows x ncols matrix A by the ncols x mcols matrix   //
//     B to form the nrows x mcols matrix C, i.e. C = A B.                    //
//     The matrix C should be declared as double C[nrows][mcols] in the       //
//     calling routine.  The memory allocated to C should not include any     //
//     memory allocated to A or B.                                            //
//                                                                            //
//  Arguments:                                                                //
//     double *C    Pointer to the first element of the matrix C.             //
//     double *A    Pointer to the first element of the matrix A.             //
//     int    nrows The number of rows of the matrices A and C.               //
//     int    ncols The number of columns of the matrices A and the           //
//                   number of rows of the matrix B.                          //
//     double *B    Pointer to the first element of the matrix B.             //
//     int    mcols The number of columns of the matrices B and C.            //
//                                                                            //
//  Return Values:                                                            //
//     void                                                                   //
//                                                                            //
//  Example:                                                                  //
//     #define N                                                              //
//     #define M                                                              //
//     #define NB                                                             //
//     double A[M][N],  B[N][NB], C[M][NB];                                   //
//                                                                            //
//     (your code to initialize the matrices A and B)                         //
//                                                                            //
//     Multiply_Matrices(&C[0][0], &A[0][0], M, N, &B[0][0], NB);             //
//     printf("The matrix C is \n"); ...                                      //
////////////////////////////////////////////////////////////////////////////////
void Multiply_Matrices(double *C, double *A, int nrows, int ncols,
                                                          double *B, int mcols) 
{
   double *pB;
   double *p_B;
   int i,j,k;

   for (i = 0; i < nrows; A += ncols, i++) 
      for (p_B = B, j = 0; j < mcols; C++, p_B++, j++) {
         pB = p_B;
         *C = 0.0; 
         for (k = 0; k < ncols; pB += mcols, k++) 
            *C += *(A+k) * *pB;
      }
}

I don't know why, I must have missed exactly this link to see through - I think it will be really helpful.
Finally it proves again: I certain code says more than a thousand links

So now I revised it to the following and it almost looks like my original code (not sure about some obfuscations, e.g. *pB and *p_B)^^

Code: Select all

void MatrixMatrixMult (double *A,  double *B, double *C, int N, int M, int K) // A[N][M] * B[M][K] = C[N][K]
{
   double *pB;
   double *p_B;
   int i,j,s;

   for (i = 0; i < N; A += M, i++)
      for (p_B = B, j = 0; j < K; C++, p_B++, j++) {
         pB = p_B;
         *C = 0.0;
         for (s = 0; s < M; pB += K, s++)
            *C += *(A+s) * *pB;
      }
}

Thanks again, I will now see how to get through the rest of my matrix algebraics

[HaWe]

works! :)

Code: Select all

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/ioctl.h>

#include "lms2012.h"
#include "ev3_button.h"
#include "ev3_lcd.h"
#include "ev3_constants.h"
#include "ev3_command.h"
#include "ev3_timer.h"
#include "ev3_lcd.h"
#include "ev3_sound.h"
#include "ev3_output.h"



//     Multiply_Matrices(&A[0][0], &B[0][0], &C[0][0], N,M,K);             //

void MatrixMatrixMult (double *A, double *B, double *C, int N, int M, int K) // A[N][M] * B[M][K] = C[N][K]
{
   double *pB;
   double *p_B;
   int i,j,s;

   for (i = 0; i < N; A += M, i++)
      for (p_B = B, j = 0; j < K; C++, p_B++, j++) {
         pB = p_B;
         *C = 0.0;
         for (s = 0; s < M; pB += K, s++)
            *C += *(A+s) * *pB;
      }
}


int main() {
  char buf[120];

  double A[2][2], B[2][2], C[2][2];
  double O[3][3], P[3][3], Q[3][3];
  
  ButtonLedInit();
  LcdInit();
  LcdClean();

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

    B[0][0]=10;  B[0][1]=30;
    B[1][0]=20;  B[1][1]=40;
    
  MatrixMatrixMult(&A[0][0], &B[0][0], &C[0][0], 2,2,2);
  sprintf (buf, "%4.0f %4.0f", A[0][0], A[0][1]); LcdText(1, 0,10, buf);
  sprintf (buf, "%4.0f %4.0f", A[1][0], A[1][1]); LcdText(1, 0,20, buf);
  LcdText(1, 0,30, "*");
  sprintf (buf, "%4.0f %4.0f", B[0][0], B[0][1]); LcdText(1, 0,40, buf);
  sprintf (buf, "%4.0f %4.0f", B[1][0], B[1][1]); LcdText(1, 0,50, buf);
  LcdText(1, 0,60, "=");
  sprintf (buf, "%4.0f %4.0f", C[0][0], C[0][1]); LcdText(1, 0,70, buf);
  sprintf (buf, "%4.0f %4.0f", C[1][0], C[1][1]); LcdText(1, 0,80, buf);


  while(ButtonWaitForAnyPress(100) != BUTTON_ID_ESCAPE)

  LcdExit();
  ButtonLedExit();
  return 1;

}

[HaWe]

...
I'm still stuck to this **A thing.
I can't imagine how to keep the definitions for the matrices like
double A[2][2], B[2][2], C[2][2];

and how to change the following (correct working) recursive determinate function to work with this definition...?

Code: Select all

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/ioctl.h>

#include "lms2012.h"
#include "ev3_button.h"
#include "ev3_lcd.h"
#include "ev3_constants.h"
#include "ev3_command.h"
#include "ev3_timer.h"
#include "ev3_lcd.h"
#include "ev3_sound.h"
#include "ev3_output.h"


// Recursive definition of determinate using expansion by minors.

double MatrixDet(double **a, int n)
{
   int i,j,j1,j2;
   double det = 0;
   double **m = NULL;

   if (n < 1) { /* Error */

   } else if (n == 1) { /* Shouldn't get used */
      det = a[0][0];
   } else if (n == 2) {
      det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
   } else {
      det = 0;
      for (j1=0;j1<n;j1++) {
         m = malloc((n-1)*sizeof(double *));
         for (i=0;i<n-1;i++)
            m[i] = malloc((n-1)*sizeof(double));
         for (i=1;i<n;i++) {
            j2 = 0;
            for (j=0;j<n;j++) {
               if (j == j1)
                  continue;
               m[i-1][j2] = a[i][j];
               j2++;
            }
         }
         det += pow(-1.0,1.0+j1+1.0) * a[0][j1] * MatrixDet(m,n-1);
         for (i=0;i<n-1;i++)
            free(m[i]);
         free(m);
      }
   }
   return(det);
}

int main() {
  char buf[120];

  double A[2][2], B[2][2], C[2][2];
  double P[3][3], Q[3][3], R[3][3];
  
  double det;
  
  ButtonLedInit();
  LcdInit();
  LcdClean();

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

    B[0][0]=10;  B[0][1]=30;
    B[1][0]=20;  B[1][1]=40;
    
  det=MatrixDet(A, 2);

  sprintf (buf, "%4.0f %4.0f", A[0][0], A[0][1]); LcdText(1, 0,10, buf);
  sprintf (buf, "%4.0f %4.0f", A[1][0], A[1][1]); LcdText(1, 0,20, buf);
  LcdText(1, 0,30, "Determinate=");

  sprintf (buf, "%4.0f ", det); LcdText(1, 0,40, buf);

  while(ButtonWaitForAnyPress(100) != BUTTON_ID_ESCAPE)

  LcdExit();
  ButtonLedExit();
  return 1;

}

[gloomyandy]

Doc,
does the code that you have posted even compile? I'd be surprised if it does. You really need to understand exactly what a declaration like:
double A[10][5]
means and how it is laid out in memory etc. You simply can't pass an array declared as above to a function that takes a parameter like double **a. It simply won't work. A 2D array of doubles in C is not an array of pointers to an array of doubles (which is what your code seems to expect it to be). It is simply a lump of memory large enough to hold (in this case) 10*5 doubles. Did you read the tutorial I posted way back at the start of this thread?

[HaWe]

Andy,
of course the code doesn't compile, there's an error about wrong pointer type - of course !
I know that the function declaration has to be changed, but I have no idea into what -
I do understand simple pointer operations but not so complicated ones.
For the moment I don't need to understand how it will work, I just need it to work, maybe in few years of practice I'll be ready to understand it - but it has to work already now.

A workaround might be to pass the arrays[][] like I did before in the multiplication function but then assign the passed-through matrices to datatypes which match the existing recursive function and then proceed to work with these transformed ones.
Or, alternatively, a different code which matches my matrix definitions.
But I have no idea how to do this, that's why Im looking for already existing tested library functions for this.

[HaWe]

Andy,
I thought about what you wrote about plain C and matrix operations.
Fact is that almost all of my source code for any of my projects is using matrices - either my Astar or my Neural Nets or my map-and explorere robots or sensor fusioning for movement vectors and Kalman filtering or my 5DOF robotic arm control.
I always have to do things like

Code: Select all

for (i=0;i<n-1;i++)       
         for (i=1;i<n;i++) {
            j2 = 0;
            for (j=0;j<n;j++) {
               if (j == j1)
                  continue;
               m[i-1][j2] = a[i][j];
               j2++;
            }
         }
         det += pow(-1.0,1.0+j1+1.0) * a[0][j1] * MatrixDet(m,n-1);

by very generally kept calculation functions at a more abstract level which are called for certain requirements by other functions for more concrete calculation levels.
To have pointers is fine just in order to point to already existing array addresses in the memory and then to manipulate the array cell values directly by the called function,
but it is even more significant to address and count through cell-by-cell and calculate one by the other as I showed above.

So if I can't pass easily arrays like A[N][M][K] to functions (by value or by reference) and then easily access each "cell" A[n][m][k] and its neighboring cells like, e.g. A[n][m-1][k+1] of a more-dim array like I'm doing it in the code above then indeed plain C is not the right programming language for these projects. I was thinking about it very closely last evening until late night and tried to transform my NXC matrix code to plain C - and failed.

I think I have to throw in the towel about plain C. Thanks for your advice.

[HaWe]

It did not leave me in peace...

after I have to pass the matrix sizes anyway (int N, int M, int K, int order) then I'm able to pass the matrices by their actual sizes as well (A[N][M], B[M][K],...)

and by sort of this it works - and it's quite easy anyway!
(in case the number of rows of the current square matrix can be calculated by int u=(int)sqrt(sizeof(R)/sizeof(R[0][0])); )

matrix multiplication:

Code: Select all

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/ioctl.h>

#include "lms2012.h"
#include "ev3_button.h"
#include "ev3_lcd.h"
#include "ev3_constants.h"
#include "ev3_command.h"
#include "ev3_timer.h"
#include "ev3_lcd.h"
#include "ev3_sound.h"
#include "ev3_output.h"





void MatrixMatrixMult(int N, int M, int K, double A[N][M], double B[M][K], double C[N][K]){
  int n, k, m;
  for (n=0; n<N; ++n) {
    for (k=0; k<K; ++k) {
      C[n][k]=0;
      for (m=0; m<M; ++m) {
        C[n][k]=C[n][k] + A[n][m]*B[m][k];
      }
    }
  }
}


int main() {
  char buf[120];

  double A[2][2], B[2][2], C[2][2];
  double P[3][2], Q[2][3], R[3][3];

  ButtonLedInit();
  LcdInit();
  LcdClean();

    P[0][0]=1;   P[0][1]=4;
    P[1][0]=2;   P[1][1]=5;
    P[2][0]=3;   P[2][1]=6;

    Q[0][0]=10;  Q[0][1]=30; Q[0][2]=50;
    Q[1][0]=20;  Q[1][1]=40; Q[1][2]=60;

  MatrixMatrixMult(3,2,3, P, Q, R);

  sprintf (buf, "%4.0f %4.0f %4.0f", R[0][0], R[0][1], R[0][2]); LcdText(1, 0,70, buf);
  sprintf (buf, "%4.0f %4.0f %4.0f", R[1][0], R[1][1], R[1][2]); LcdText(1, 0,80, buf);
  sprintf (buf, "%4.0f %4.0f %4.0f", R[2][0], R[2][1], R[2][2]); LcdText(1, 0,90, buf);

  while(ButtonWaitForAnyPress(100) != BUTTON_ID_ESCAPE)

  LcdExit();
  ButtonLedExit();
  return 1;

}

matrix determinant:

Code: Select all

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <string.h>
#include <sys/mman.h>
#include <sys/ioctl.h>

#include "lms2012.h"
#include "ev3_button.h"
#include "ev3_lcd.h"
#include "ev3_constants.h"
#include "ev3_command.h"
#include "ev3_timer.h"
#include "ev3_lcd.h"
#include "ev3_sound.h"
#include "ev3_output.h"


double MatrixDet(int N, double A[N][N])
{
    int i,j,i_count,j_count, count=0;
    double Asub[N-1][N-1], det=0;
    
    if(N==1) return A[0][0];
    if(N==2) return (A[0][0]*A[1][1] - A[0][1]*A[1][0]);

    for(count=0; count<N; count++)
    {
        i_count=0;
        for(i=1; i<N; i++)
        {
            j_count=0;
            for(j=0; j<N; j++)
            {
                if(j == count) continue;
                Asub[i_count][j_count] = A[i][j];
                j_count++;
            }
            i_count++;
        }
        det += pow(-1, count) * A[0][count] * MatrixDet(N-1,Asub);
    }
    return det;
}


int main() {

  ButtonLedInit();
  LcdInit();
  LcdClean();
  
  char   buf[120];
  double det;

  double A[2][2], B[2][2], C[2][2];
  double P[3][2], Q[2][3], R[3][3];

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

    B[0][0]=10;  B[0][1]=30;
    B[1][0]=20;  B[1][1]=40;
    
    P[0][0]=1;   P[0][1]=4;
    P[1][0]=2;   P[1][1]=5;
    P[2][0]=3;   P[2][1]=6;

    Q[0][0]=10;  Q[0][1]=30; Q[0][2]=50;
    Q[1][0]=20;  Q[1][1]=40; Q[1][2]=60;
    
    R[0][0]=11;  R[0][1]=44; R[0][2]=77;
    R[1][0]=22;  R[1][1]=55; R[1][2]=88;
    R[2][0]=34;  R[2][1]=67; R[2][2]=99;

  det=MatrixDet(2, A);
  sprintf (buf, "det A 2x2=%4.0f", det); LcdText(1, 0,60, buf);
  
  det=MatrixDet(3, R);
  sprintf (buf, "det R 3x3=%4.0f", det); LcdText(1, 0,80, buf);

  while(ButtonWaitForAnyPress(100) != BUTTON_ID_ESCAPE)

  LcdExit();
  ButtonLedExit();
  return 1;

}

so as this matrix thing is fixed....

...now what about daisy-chaining?