Help with cblas_dgemm matrix multiplication

[SpiffyEh]

I'm having problems running cblas_dgemm on a matix matrix multiplication.
I have the following matricies

double * mass = new double[n];                  
double (* pos)[NDIM] = new double[n][NDIM];
double tempPos[NDIM];
double tempMass[nlocal];
double mass_avg[1];
double pos_avg[NDIM];

The idea is to multiply the mass * tempMass and put it in mass_avg
and multiple pos * tempPos and put it in pos_avg
nlocal is less than n and it needs to multiple up to nlocal rows not n

I'm using the following call for mass and it works as expected:

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, 1, nlocal, 1.0, mass, nlocal, tempMass, 1, 1.0, m_ave, 1);

It's the position call that isn't working for me. I have the following for that one:

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nlocal, 1, NDIM, 1.0, pos[0], NDIM, tempPos, 1, 0.0, pos_ave, 1);

I've been struggling with this for hours and can't figure out why the position call isn't working. I keep getting a segmentation fault. Any ideas about what's going on? Or what I'm doing wrong in this situation?

 

[Mark44]

I'm having problems running cblas_dgemm on a matix matrix multiplication.
I have the following matricies

double * mass = new double[n];                  
double (* pos)[NDIM] = new double[n][NDIM];
double tempPos[NDIM];
double tempMass[nlocal];
double mass_avg[1];
double pos_avg[NDIM];

The idea is to multiply the mass * tempMass and put it in mass_avg
and multiple pos * tempPos and put it in pos_avg
nlocal is less than n and it needs to multiple up to nlocal rows not n

I'm using the following call for mass and it works as expected:

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, 1, nlocal, 1.0, mass, nlocal, tempMass, 1, 1.0, m_ave, 1);

It's the position call that isn't working for me. I have the following for that one:

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nlocal, 1, NDIM, 1.0, pos[0], NDIM, tempPos, 1, 0.0, pos_ave, 1);

Here is some documentation for cblas_dgemm that I found: (http://www.psatellite.com/matrixlib/api/lapack.html ). Make sure that what you are passing in each parameter of this function matches exactly with the description in the table. You don't show some of the values in the call that's failing, so I can't check that the values are OK.

Do you understand what each parameter represents?

I've been struggling with this for hours and can't figure out why the position call isn't working. I keep getting a segmentation fault. Any ideas about what's going on? Or what I'm doing wrong in this situation?

 

[SpiffyEh]

I understand what the parameters mean. I just don't see what I'm doing wrong. I have checked all of the meanings multiple times and the values I'm passing.

the mass matrix is n (303) x 1
pos is n(303) x NDIM (3)
tempPos is NDIM x 1
tempMass is nlocal(300) x 1
mass_avg is 1x1
pos_avg is NDIM(3) x 1

n is actually 3 more than nlocal. I included all of the values up above.
I don't think I'm passing any of the values in wrong which is why I'm so confused.

 

[Mark44]

I understand what the parameters mean. I just don't see what I'm doing wrong. I have checked all of the meanings multiple times and the values I'm passing.

the mass matrix is n (303) x 1

No, it is 1 X 303, meaning that it has 1 row, with 303 columns per row. Strictly speaking mass is a one-dimension array with 303 elements.

The same goes for all of your other one-dimension arrays, including tempPos, tempMass, and pos_avg. I don't know why you made mass_avg an array, since it has only one element, but what you did seems to have worked for you. An alternate way to do things would be to declare mass_ave as a double, and then just pass &mass_ave as a parameter.

pos is n(303) x NDIM (3)
tempPos is NDIM x 1
tempMass is nlocal(300) x 1
mass_avg is 1x1
pos_avg is NDIM(3) x 1

pos is actually 3 x 303, based on what I think you wanted it to be. Your declaration is fouled up, and I think this is what is causing your problem. See below.

n is actually 3 more than nlocal. I included all of the values up above.
I don't think I'm passing any of the values in wrong which is why I'm so confused.

I believe the problem is this line:

double (* pos)[NDIM] = new double[n][NDIM];

Your declaration for pos is not right. On the left side, you are telling the compiler that pos is a pointer to an array with 3 elements. On the right side, you are allocating space on the heap for an array of 303 elements, each element of which is a 3 element array. In other words, the right side allocates space for a 303 x 3 matrix rather than a 3 x 303 matrix.

If you change the declaration above to the following, I think that will fix your problem.

double (* pos)[NDIM] = new double[NDIM][n];

Alternatively, you could declare pos this way, allocating space on the stack:

double pos[NDIM][n];

This would give you an uninitialized matrix with 3 rows and 303 columns. pos[0] would evaluate to the address of the first row (which is also the address of pos[0][0]).

 

[SpiffyEh]

all of the arrays i said the dimensions for are in row major or are supposed to be. I can't change the pos array dimensions or the rest of the program will break. It's written by someone else as a simulation and is pretty long, I'm just adding averages and a couple other things to it.

I can change anything besides the mass and pos arrays. The position matrix should be 303 rows by 3 columns. Am I doing the call to dgemm right in this situation?

I didn't realize I could pass in a double for mass in that way, thanks for letting me know!

 

[Mark44]

If pos is a matrix with 303 rows, and 3 columns per row, then it should be declared in either of these ways:

double pos[n][NDIM]; // Uninitialized matrix allocated on stack.
OR
double (*pos)[n] = new double[n][NDIM]; // Uninitialized matrix allocated on heap.

What you showed in post #1 is wrong..

double (* pos)[NDIM] = new double[n][NDIM];

What I'm calling a "matrix" above is really a two-dimension array of arrays. You can think of pos as being an array of 303 elements (the rows), where each element(or row) is an array with three elements (the columns).


From what you have already said, you want to do this multiplication -- pos * tempPos -- and store the result in pos_avg.

pos is 303 x 3, tempPos is 3 x 1, so the product will be 303 x 1, not 3 x 1 as you have for pos_avg. You need to declar tempPos to be the right size: double pos_avg[303].

I also think that a couple of the parameters in your call to cblas_dgemm are not right.

Parameters 1, 2, and 3 are OK, I think.

The following applies to parameters 4 through 12.
4. nlocal = #rows in A matrix == this should be n (303), not nlocal (300)
5. 1 = #cols in B == OK
6. NDIM = #cols in A == OK
7. 1.0 - a scalar == OK, I guess. I don't know what this is for.
8. pos[0]The array for A == OK, once you get your array declared right.
9. NDIM = stride for A == OK. The stride is the number of elements per row.
10.tempPos = The array for B == OK
11. 1 = stride for B == OK
12. 0.0 - scalar to multiply by == ?? Seems like it should be 1.0
13. pos_ave = the array for C == OK
14. 1 = stride for C == OK

 

[SpiffyEh]

Thanks for the help, I appreciate you spending time on trying to figure this out. I ended up figuring out that the pos matrix needed to be translated so the second parameter was CblasTrans and I needed to make tempPos larger and it worked.

 

[Mark44]

Glad I could help out.