Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Trouble with BLAS subroutine dgemv ->vector- matrix product

  1. May 15, 2010 #1
    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>

    // [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

  2. jcsd
  3. May 15, 2010 #2
    Re: Trouble with BLAS subroutine dgemv -->vector- matrix product

    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)
  4. May 15, 2010 #3
    Re: Trouble with BLAS subroutine dgemv -->vector- matrix product

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

    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,


    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.
  5. May 15, 2010 #4
    Re: Trouble with BLAS subroutine dgemv -->vector- matrix product


    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
  6. May 15, 2010 #5


    Staff: Mentor

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

    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.
  7. May 15, 2010 #6
    Re: Trouble with BLAS subroutine dgemv -->vector- matrix product


    I just tried it with radians, and it works. mystery solved.
  8. May 15, 2010 #7
    Re: Trouble with BLAS subroutine dgemv -->vector- matrix product


    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?


    -miss fangula
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Similar Discussions: Trouble with BLAS subroutine dgemv ->vector- matrix product