# Thread: matrix mult - double ptr in col major

1. ## 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. >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.

double **matrixA, **matrixB, **matrixC;

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

3. i am not sure what you are trying to do but you can use mutli-dimentional arrays to store whole tables.

4. 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
*/```

5. After reading Dave's link, ignore my post.

6. 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]+ Y = 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. 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...
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;
}```

8. 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. 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.
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.

10. 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. 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
*/```

12. 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. 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.

14. 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