Trouble with BLAS subroutine dgemv ->vector- matrix product

[missfangula]

Trouble with BLAS subroutine dgemv -->vector- matrix product

Hello all,

So, I am trying to write some code that will compute the vector-matrix product using the BLAS subroutine dgemv. I have written the code below, but I am getting really wacky numbers! You will see that it takes array "Array" as the known matrix , and vector 'vector' as known vector, and outputs a product vector.

Does anyone know what maybe askew?

Thanks a million,
miss fangula


#import <Foundation/Foundation.h>
#include <Accelerate/Accelerate.h>
#include <stdlib.h>
#include <stdio.h>
//#include <math.h>

//*******FINDS X IN MATRIX EQUATION AB = X, WHERE A AND B ARE KNOWN MATRICES*******//
//
// [a a a] [b b b] [x x x]
// [a a a]*[b b b] = [x x x] {x,x,x,x,x,x,x,x,x} = ?
// [a a a] [b b b] [x x x]
//
//*********************************************************************************//

void printMatrix(char* name, int n, double* b, int lda);

/*Define all subroutine parameters*************************/
#define M 4 //Define 'Array' matrix size and number of rows in 'vector'
#define N 4 //Define number of columns in 'Array' (if a vector, then probably 1)
#define ALPHA 4 //Scaling constant alpha
#define LDA 4 //Leading dimension of 'Array'
#define INCX 1 //Stride value for 'vector' and can have any value
#define BETA 0 //Scaling constant beta
#define INCY 1 //Stride for 'y'
/**********************************************************/

int main()
{
double *vector;
double *y;
int m = M;//Define local variable n to be N
int n = N;//
int alpha = ALPHA;
int lda = LDA; //Define local variable lda to be LDA
int incx = INCX;
int beta = BETA;
int incy = INCY;

//**********
double a = 0.0;
double b = 180.0;
double R = 1.54;
//**********


double Array[16] = //Declare a NxN array 'Array'
{
-cos(a), sin(a)*cos(b), sin(a)*sin(b), 0.0,
-sin(a), -cos(a)*cos(b), -cos(a)*sin(b), 0.0,
0.0, -sin(b), cos(b), 0.0,
-R*cos(a), R*sin(a)*cos(b), R*sin(a)*sin(b), 1.0
};




vector = calloc(4, sizeof (double)); //Declare a vector to as other unknown in equation (see algorithm above)
y = calloc(4, sizeof (double));

*vector = 0.0; //Define first element of vector b
*(vector + 1) = 0.0; //Define second element of vector b
*(vector + 2) = 0.0;
*(vector + 3) = 1.0;


cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, Array, lda,
vector, incx, beta, y, incy);
/*cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1.0, m, 3,
x, 1, 0.0, y, 1);*/

printMatrix("Answer is: ", n, y, lda); //Print the solution matrix 'x'

free(vector); //free the input vector
free(y); //Free the solution vector

return 666; //since int main function, return some integer
}

void printMatrix (char* name, int n, double* b, int lda) //Function to print the solution matrix
{
int loopNumber;
int numTimes = n;
printf("\n %s\n", name); //Print name of solution matrix 'x'
for (loopNumber = 0; loopNumber < numTimes; loopNumber++) //For each element of vector
{
printf("%6.2f", b[loopNumber]); //Print the respective element at loopNumber position
printf( "\n" ); //Space
}

}

 

[missfangula]

i should add that i am running this on a mac. so i am using the accelerate framework for lapack (and blas, too, I believe)

 

[TheoMcCloskey]

Check your use of sin(b) and cos(b). Those functions are expecting radian arguments.

Also,
What numbers were you expecting? What numbers are you getting?

Be aware, the routine "dgemv" returns a matrix-vector operation

y := alpha*A*x + beta*y,

or

y := alpha*A'*x + beta*y,


depending on the value of "CblasNoTrans". Here, "alpha" and "beta" are scalar.


I suggest you try a simple test case with some know variables.

 

[missfangula]

Hi,

I set beta scalar to be 0.0, and alpha scalar to be 1.0, so now it should output the product of a matrix A and a vector x, correct? A*x = ?

I tested the program with various 2x2, 3x3, and 4x4 matrices, all with constants as elements in the array. I got the correct answers.

Then I tried it as posted in this forum, with variables for elements in the array. I get wrong answers. I checked the results with an online matrix multiplication calculator (joemath.com)

I think I have narrowed down the problem to either one of two things:
1. I cannot use variables for the array elements
2. sine is in radians, and I am putting in degrees

I will now try out the radian version, and see what happens. Otherwise, does anyone know of a rule that you cannot use variables as array elements in c? I thought that was only when allocating the size of the array (need a constant for array size)...

-miss fangula

 

[Mark44]

It shouldn't be a problem using variables to intitialize the array, but it's definitely a problem that you are using degrees instead of radians.

 

[missfangula]

hi,

I just tried it with radians, and it works. mystery solved.

 

[missfangula]

Hi

I had another question, kind of silly maybe. What is the difference between a matrix-vector product and a matrix-matrix product in programming? Aren't they are just a bunch of arrays? If s0, then can I use the dgemv subroutine to multiply two 4x4 matrices, do I have to use another subroutine, like dgemms?

Thanks!

-miss fangula