Thread: matrix mult - double ptr in col major

  1. #1
    Registered User
    Join Date
    Mar 2005
    Posts
    17

    matrix mult - double ptr in col major

    Hey,

    I'm creating a matrix-matrix multiplier in C that will use the same format as the GOTO BLAS and ATLAS implementations. I need to allocate memory and store the matrices in COL MAJOR format not in row major as is the standard C representation.

    Can you tell me if I am on the right track here:
    Code:
    //global vars
    double **matrixA, **matrixB, **matrixC;
    
    int main(int argc, char *argv[])
    {
    	int i, c, row, col, j;
    	int n_rows = 10; //ASSUME SQUARE FOR NOW!
    	int m_cols =  10;
    	int mbyn = m_cols * n_rows;
            int i, j;
            double n = 0;
    
    //allocate the number of cols
    matrixA = (double**)malloc( sizeof( double* ) * m_cols); 
    
    //allocate the overall size??
    matrixA[0] = (double *) malloc( sizeof( double*) * mbyn);
    	
    for (i=1; i<m_cols; i++)	
    {       //for each col allocate the rows
    	matrixA[i] = *matrixA + i * n_rows; 
    }
    I then have a function fill matrix so:

    Code:
      for (i=0; i<SIZE; i++) //
        for (j=0; j<SIZE; j++)
          matrix[j][i] = n++;  //or a[i][j];  //fill it col,row - like a graphics point

    I think the malloc is correct but i'm not so sure about the filling of the matrix. Is this right??

    I want it to work like this:
    Image Attachment
    (sorry don't know how to display the image in here.)

    So that I have one long linear array, that stores one row, then the next row etc. This means that I can skip to the next element one by one the whole way throught the matrix. To get to the next row of a matrix with n cols I would simply move n elements from the start pointer.

    thanks for your help,

    Colm Mitchell

  2. #2
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    >So that I have one long linear array, that stores one row, then the next row etc.

    The only way to do that is just have a single pointer, and access column and row using pointer arithmetic.

    Instead of:
    double **matrixA, **matrixB, **matrixC;

    You would have:
    double *matrixA, *matrixB, *matrixC;

  3. #3
    Registered User mrafcho001's Avatar
    Join Date
    Jan 2005
    Posts
    483
    i am not sure what you are trying to do but you can use mutli-dimentional arrays to store whole tables.

  4. #4
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    Second example here.

    I know this came up here recently and Salem had a really nice solution... this thread: here.

    [edit=3]I'm not sure if this is quite what you're after.
    Code:
    #include <stdio.h>
    #include <stdlib.h>
    
    int **foo(int nrows, int ncolumns)
    {
       int i, **array2 = malloc(nrows * sizeof *array2);
       if ( array2 )
       {
          array2[0] = malloc(nrows * ncolumns * sizeof *array2[0]);
          if ( array2[0] )
          {
             for ( i = 1; i < nrows; i++ )
             {
                array2[i] = array2[0] + i * ncolumns;
             }
          }
       }
       return array2;
    }
    
    void bar(int **array)
    {
       free(array[0]);
       free(array);
    }
    
    int main(void)
    {
       int r, c, i = 0, rows = 3, cols = 5, **myarray = foo(cols, rows);
       for ( c = 0; c < cols; ++c )
       {
          for ( r = 0; r < rows; ++r )
          {
             myarray[c][r] = i++;
             printf("%p : myarray[%d][%d] = %d\n",
                    (void*)&myarray[c][r], c, r, myarray[c][r]);
          }
       }
       bar(myarray);
       return 0;
    }
    
    /* my output
    007A3208 : myarray[0][0] = 0
    007A320C : myarray[0][1] = 1
    007A3210 : myarray[0][2] = 2
    007A3214 : myarray[1][0] = 3
    007A3218 : myarray[1][1] = 4
    007A321C : myarray[1][2] = 5
    007A3220 : myarray[2][0] = 6
    007A3224 : myarray[2][1] = 7
    007A3228 : myarray[2][2] = 8
    007A322C : myarray[3][0] = 9
    007A3230 : myarray[3][1] = 10
    007A3234 : myarray[3][2] = 11
    007A3238 : myarray[4][0] = 12
    007A323C : myarray[4][1] = 13
    007A3240 : myarray[4][2] = 14
    */
    Last edited by Dave_Sinkula; 03-21-2005 at 03:36 PM. Reason: [1] Found thread. [2] Isolated post. [3] Added code.
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  5. #5
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    After reading Dave's link, ignore my post.

  6. #6
    Registered User
    Join Date
    Mar 2005
    Posts
    17
    Ok - i got a bit confused by that so had to ask my supervisor again...

    Swoopy you were on the right track - I need to have a single pointer instead of a double so
    You would have:
    double *matrixA, *matrixB, *matrixC;
    would work.

    This means that when i go to fill the matrix for example i need to do it not like
    PHP Code:
    matrix[1][2] = X
    but as
    PHP Code:
    matrix[0]+ X
    where Y is the Yth element of that row. Meaning that I start a particular column and then move along the array to find the element in the Yth row of that col.

    Dave in what way will I need to change the memory allocation because of this?

    Thanks for your help guys!

    Colly.

  7. #7
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    Quote Originally Posted by collymitch
    Dave in what way will I need to change the memory allocation because of this?
    Then look at the third example there...
    Quote Originally Posted by Dave_Sinkula
    Second example here.
    To do it column major, you'd do the indexing differently. Maybe something like this.
    Code:
    #include <stdio.h>
    #include <stdlib.h>
    
    int main(void)
    {
       size_t rows = 3, cols = 5;
       double *myarray = malloc(rows * cols * sizeof *myarray);
       if ( myarray )
       {
          size_t r, c;
          for ( r = 0; r < rows; ++r )
          {
             for ( c = 0; c < cols; ++c )
             {
                size_t i = c * rows + r;
                myarray[i] = (r + 1) + (c + 1) / 10.0;
                printf("{%d,%d} %p : myarray[%2d] = %g\n", (int)r, (int)c,
                       (void*)&myarray[i], (int)i, myarray[i]);
             }
          }
          free(myarray);
       }
       return 0;
    }
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  8. #8
    Registered User
    Join Date
    Mar 2005
    Posts
    17
    Ok Dave, Thanks, that seems to work fine.

    I was unsure how you got this part:

    Code:
    size_t i = c * rows + r;
    I went thorugh it step by step and i see what it's doing but i'm not completely sure why that works. Thing is Ii think I need to know it so i can start to multiply the matrices.

    I'm going to be doing the same sort of loops aren't I?

    I had got the standard loops:
    Code:
      for (i=0; i<SIZE; i++)
        for (j=0; j<SIZE; j++)
    	{
           matrixC[i][j]=0;
           for (k=0; k<SIZE; k++)
     		printf("\n\n");
    		matrixC[i][j] += matrixA[i][k]*matrixB[k][j];
    but now these are no use are they??

    Do i need to be doing the same calculation to find i, then a similar calculation to find j??

    Thanks again,

    Clearly lost Colly.

  9. #9
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    For my own refresher, let's see if my understanding is correct.

    You want a "two-dimensional" array, but you are intending to use a library that receives a simulated 2D array that is actually a 1D array.

    This was mentioned in the c.l.c. FAQ:
    If the double indirection implied by the above schemes is for some reason unacceptable, you can simulate a two-dimensional array with a single, dynamically-allocated one-dimensional array:
    Code:
    int *array3 = (int *)malloc(nrows * ncolumns * sizeof(int));
    However, you must now perform subscript calculations manually, accessing the i,jth element with array3[i * ncolumns + j].
    But the subscripting is shown for the "normal" row-major version familiar to C. And here I believe the i and j are row and column, respectively. Since your library expects data column-major, I changed this to:
    Code:
    c * rows + r
    This would have been array3[j * nrows + i] from the c.l.c. FAQ.
    Quote Originally Posted by collymitch
    but now these are no use are they??

    Do i need to be doing the same calculation to find i, then a similar calculation to find j??
    It is just the [i][j] subscripting that needs to be changed to the [j * nrows + i] subscripting. And for [i][k] and [k][j].

    But bear with me, I find column-major confusing.
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  10. #10
    Registered User
    Join Date
    Mar 2005
    Posts
    17

    Thumbs up

    Yes that's right - i'm comparing my result with the Blas implementation which uses the DGEMM routine as described Here.

    Ah-ha. That makes sense to me now.

    So i'm gonna have

    Code:
    void mmult(void)
    {
      int i,j,k;
      int row, col;
    //#pragma omp for   
      for ( row = 0; row < rows; ++row )  //for (i=0; i<SIZE; i++)
    	for ( col = 0; col < cols; ++col )  //  for (j=0; j<SIZE; j++)
    	{	
    		i = col * rows + row;
    		matrixC[i] = 0;         //matrixC[i][j]=0;
    		for (k=0; k<rows; k++)
    		{
    			matrixC[i] += 
                                matrixA[k * rows + row] * matrixB[col * rows + k];
    			//matrixC[i][j] += matrixA[i][k]*matrixB[k][j];
    		} //end for
    	}// end for
    }
    I think that's actually working, and in Column major! > it may need to be changed to accomodate non-square matrices - i think i only need to change the ROWS in the K loop to the M value?

    You find col-major confusing - how do you think I feel haha. My supervisor told me to use f2c on a fortran program and figure it out from that!

    Thanks again Dave, you've been a great help.

    Colly.

  11. #11
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    Quote Originally Posted by collymitch
    it may need to be changed to accomodate non-square matrices - i think i only need to change the ROWS in the K loop to the M value?
    Like I said, I get a little turned around with it; here's a row-major thingy to compare to.
    Code:
    #include <stdio.h>
    #include <stdlib.h>
    
    /*
     * a[m][n] x b[n][p] = c[m][p]
     */
    void mult(const double *a, const double *b, double *c, size_t m, size_t n, size_t p)
    {
       size_t i, j, k;
       for ( i = 0; i < m; ++i )
       {
          for ( j = 0; j < p; ++j )
          {
             c [ i * p + j ] = 0;
             for ( k = 0; k < n; ++k )
             {
                c [ i * p + j ] += a [ i * n + k ] * b [ k * p + j ];
             }
          }
       }
    }
    
    void show(const double *array, size_t rows, size_t cols)
    {
       size_t r, c;
       for ( r = 0; r < rows; ++r )
       {
          for ( c = 0; c < cols; ++c )
          {
             printf("%4g ", array[r * cols + c]);
          }
          putchar('\n');
       }
       putchar('\n');
    }
    
    int main(void)
    {
       double a2x4[] = {1,2,3,4,5,6,7,8};
       double b4x3[] = {1,2,3,4,5,6,7,8,9,10,11,12};
       double c2x3[6];
       show(a2x4, 2, 4);
       show(b4x3, 4, 3);
       mult(a2x4, b4x3, c2x3, 2, 4, 3);
       show(c2x3, 2, 3);
       return 0;
    }
    
    /* my output (check http://www.mai.liu.se/~halun/matrix/matrix.html)
       1    2    3    4
       5    6    7    8
    
       1    2    3
       4    5    6
       7    8    9
      10   11   12
    
      70   80   90
     158  184  210
    */
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  12. #12
    Registered User
    Join Date
    Mar 2005
    Posts
    17
    Mmmm....

    I modified it to take in the a[m][n] x b[n][p] = c[m][p] sizes.

    I've not got it just right, it's close but close is wrong so...

    Code:
    void mmult(size_t m_rows, size_t n_common, size_t p_cols)
    {
      int i,j,k;
      int row, col;
    //#pragma omp for  
    for ( row = 0; row < m_rows; ++row ) //for (i=0; i<SIZE; i++)
    
    for ( col = 0; col < p_cols; ++col ) // for (j=0; j<SIZE; j++) {
    i = col * m_rows + row; matrixC[i] = 0; //matrixC[i][j]=0; for (k=0; k < n_common; k++) {
    matrixC[i] += matrixA[k * m_rows + row] * matrixB[col * m_rows + k];
    }
    }// end for
    }
    The red bit must not be right then. I tried to replace the row vars with the variable passed in. I bet it's something really silly.

    Output is:
    Code:
     70,  53,  80
    158, 149, 184
    Thanks,

  13. #13
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    Code:
    matrixC[i] += matrixA[k * m_rows + row] * matrixB[col * m_rows + k];
    I think m_rows should be n_common.
    Code:
    #include <stdio.h>
    
    /*
     * a[arows][acols] x b[brows][bcols] = c[arows][bcols]
     */
    void mult(const double *a, const double *b, double *c,
              size_t arows, size_t acols, size_t brows, size_t bcols)
    {
       size_t i, j, k, x, y, z;
       if ( acols != brows )
       {
          return;
       }
       for ( i = 0; i < arows; ++i )
       {
          for ( j = 0; j < bcols; ++j )
          {
             x = j * arows + i;
             c [ x ] = 0;
             for ( k = 0; k < acols; ++k )
             {
                /*
                matrixC[i] += matrixA[k * m_rows + row] * matrixB[col * m_rows + k];
                */
                y = k * arows + i;
                z = j * brows + k;
                c [ x ] += a [ y ] * b [ z ];
             }
          }
       }
    }
    
    void show(const double *array, size_t rows, size_t cols)
    {
       size_t r, c;
       for ( r = 0; r < rows; ++r )
       {
          for ( c = 0; c < cols; ++c )
          {
             size_t i = c * rows + r;
             printf("%4g ", array[i]);
          }
          putchar('\n');
       }
       putchar('\n');
    }
    
    int main(void)
    {
       double A_4Cx2R[] = {1,5,2,6,3,7,4,8};
       double B_3Cx4R[] = {1,4,7,10,2,5,8,11,3,6,9,12};
       double C_3Cx2R[6];
       show(A_4Cx2R, 2, 4);
       show(B_3Cx4R, 4, 3);
       mult(A_4Cx2R, B_3Cx4R, C_3Cx2R, 2, 4, 4, 3);
       show(C_3Cx2R, 2, 3);
       return 0;
    }
    
    /* my output (check http://www.mai.liu.se/~halun/matrix/matrix.html)
       1    2    3    4
       5    6    7    8
    
       1    2    3
       4    5    6
       7    8    9
      10   11   12
    
      70   80   90
     158  184  210
    */
    I'd thrown debugging printfs to help me get things straight.
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  14. #14
    Registered User
    Join Date
    Mar 2005
    Posts
    17

    Smile

    Excellent. That seems to be working now. I could have sworn that I'd tried that! maybe i didn't save it right or something. thanks Dave you've been a great help.

    That makes it easier to read and that's a good tip too.

    Thanks,

    Colly.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. C - access violation
    By uber in forum C Programming
    Replies: 2
    Last Post: 07-08-2009, 01:30 PM
  2. Compiling c++ intrinsics commands
    By h3ro in forum C++ Programming
    Replies: 37
    Last Post: 07-13-2008, 05:04 AM
  3. Copying 2-d arrays
    By Holtzy in forum C++ Programming
    Replies: 11
    Last Post: 03-14-2008, 03:44 PM
  4. Conversion From C++ To C
    By dicon in forum C++ Programming
    Replies: 2
    Last Post: 06-10-2007, 02:54 PM
  5. Hi, could someone help me with arrays?
    By goodn in forum C Programming
    Replies: 20
    Last Post: 10-18-2001, 09:48 AM