Matrix-vector product with dgemm/dgemv

Question

Using Lapack with C++ is giving me a small headache. I found the functions defined for fortran a bit eccentric, so I tried to make a few functions on C++ to make it easier for me to read what's going on.

Anyway, I'm not getting th matrix-vector product working as I wish. Here is a small sample of the program.

smallmatlib.cpp:

#include <cstdio>
#include <stdlib.h>


extern "C"{
    // product C= alphaA.B + betaC                                               
   void dgemm_(char* TRANSA, char* TRANSB, const int* M,
               const int* N, const int* K, double* alpha, double* A,
               const int* LDA, double* B, const int* LDB, double* beta,
               double* C, const int* LDC);
    // product Y= alphaA.X + betaY                                               
   void dgemv_(char* TRANS, const int* M, const int* N,
               double* alpha, double* A, const int* LDA, double* X,
               const int* INCX, double* beta, double* C, const int* INCY);
   } 

void initvec(double* v, int N){
    for(int i= 0; i<N; ++i){
        v[i]= 0.0;
        }
   }

void matvecprod(double* A, double* v, double* u, int N){ 
    double alpha= 1.0, beta= 0.0;
    char no= 'N', tr= 'T';
    int m= N, n= N, lda= N, incx= N, incy= N;
    double* tmp= new double[N];
    initvec(tmp, N);
    dgemv_(&no,&m,&n,&alpha,A,&lda,v,&incx,&beta,tmp,&incy);
    for(int i= 0; i<N; ++i){
        u[i]= tmp[i];
        }
    delete [] tmp;
    }

void vecmatprod(double* v, double* A, double* u, int N){
    double alpha= 1.0, beta= 0.0;
    char no= 'N', tr= 'T';
    int m= N, n= 1, k= N, lda= N, ldb= N, ldc= N;
    double* tmp= new double[N];
    initvec(tmp, N);
    dgemm_(&no,&no,&m,&n,&k,&alpha,A,&lda,v,&ldb,&beta,tmp,&ldc);
    for(int i= 0; i<N; ++i){ 
        u[i]= tmp[i];
        }
    delete [] tmp;
    }

smallmatlib.h:

#ifndef SMALLMATLIB_H
#define SMALLMATLIB_H

void initvec(double* v, int N);

void matvecprod(double* A, double* v, double* u, int N);

void vecmatprod(double* v, double* A, double* u, int N);

#endif

smallmatlab.cpp:

#include "smallmatlib.h"
#include <cstdio>
#include <stdlib.h>
#define SIZE 2

int main(){
  double A[SIZE*SIZE]=
    { 1,2,
      3,4 };
  double v[SIZE]= {2,5.2};
  double u[SIZE]= {0,0};
  matvecprod(A,v,u,SIZE);
  printf("%f %f\n",u[0],u[1]);
  vecmatprod(v,A,u,SIZE);
  printf("%f %f\n",u[0],u[1]);
  return 0;
  }

Compiling:
g++ -c smallmatlib.cpp
g++ smallmatlab.cpp smallmatlib.o -L/usr/local/lib -lclapack -lcblas

Now the function matvecprod is the problem. With the example matrix A and example vector v, it should produce an output like

12.4..    26.8..

but instead it prints out

2.00..    0.00..

I tried to produce the correct result with both dgemm and dgemv, but wasn't able to. I have a hunch that my variables incx and incy do not have correct values, but it's difficult to find an explanation for them that I'd understand.

A smaller problem is that at the moment I can't use them like vecmatprod(v,A,v,SIZE) - that is, I always have to define the vector u, that will hold the result, separately, and call vecmatprod(v,A,u,SIZE). Any help would be appreciated.

As a side note, I'm also a beginner at C++, so I appreciate any criticism/suggestion you might have about my code.

Answer 1

You are right, the problem is in incx value - it should be 1, take a look at reference.

INCX is INTEGER
  On entry, INCX specifies the increment for the elements of
  X. INCX must not be zero.

So this value should be used when values of vector x is not placed one by one (vector of complex values for example, when you want to use only real part).

Also you can not use vecmatprod(v,A,v,SIZE) with v as both x and y. This is because how matrix-vector multiplication works (see wikipedia for example). You need values of original x the whole time to produce correct result. Small example:

y = A * x 

where

y = [ y1 y2 ]
A = [ [a11 a12] [a21 a22] ]
x = [ x1 x2 ]

And we calculate y like this

y1 = a11 * x1 + a12 * x2
y2 = a21 * x1 + a22 * x2

You can see, that when we calculate y2 we need x1 and x2, but if you use x = A * x(without temporary vector) you will replace x1 with y1 thus produce wrong answer.