/*********************************************************************/
/* Multiplication banded-matrix x vector: Data structure and Blas
/*  Compile with:
/ *  gcc -O3 gbmv.c /usr/lib/libblas.so.3
/*********************************************************************/
#include <stdio.h>
#include <time.h>

#define N 80 
#define KB 2
#define NIT 2000

void dgbmv_(char*, int*, int*, int*, int*, double*, double*, int*, double*, int*, double*, double*, int*);

/* Structure pour matrice bande */
typedef struct
{  int nA;                 /* ordre de la matrice*/
   int kb;                 /* 1+2kb bandes non-nulle */
   int ib[KB+1];   
   double ab[1+2*KB][N*N];  
}  Matb;
int main()
{  
  int i,j,k;
  Matb A;  
  int n=N, nn=n*n, kb=KB, inc=1, lda=2*n+1, ku=n, nit=NIT;
  double alpha=1, beta=0;
  double x[N*N],y[N*N];
  double A1[N*N][2*N+1];
  char trans='N';
 
  clock_t clock0, clock1;

  /*Initialisation of the matrix and vectors */
  
  A.nA=nn; A.kb=kb;
  A.ib[0]=n; A.ib[1]=1; A.ib[2]=0;

  for (i=0; i<A.nA; i++)
  {  x[i]=1; y[i]=0;
     for (j=0; j<5; j++) A.ab[j][i]=0;
  }
  
  for (j=0; j<A.nA; j++) A.ab[0][j]=-1;
  for (j=0; j<A.nA; j++) A.ab[1][j]=-1;
  for (j=0; j<A.nA; j++) A.ab[2][j]=4;
  for (j=0; j<A.nA; j++) A.ab[3][j]=-1;
  for (j=0; j<A.nA; j++) A.ab[4][j]=-1;
  
/* Bandwise multiplication
**************************/
  
  clock0 = clock();

  for (k=0; k<nit; k++)
  {  for (i=0; i<3; i++)
       for (j=A.ib[i]; j<A.nA; j++) y[j-A.ib[i]]+=A.ab[i][j]*x[j];
     
     for (i=3; i<5; i++)
       for (j=0; j<A.nA-A.ib[i-3];j++) y[A.ib[i-3]+j]+=A.ab[i][j]*x[j];
  }
     
  clock1 = clock();

  if (N<6) { for (i=0; i<A.nA; i++)  printf(" %e\n",y[i]);};
  printf("\n Band multiplication CPU-time :  %e miliseconds\n",
	 (1.0e+3*(double)(clock1 - clock0))/CLOCKS_PER_SEC);

 /* Blas
*********/
 
  for (i=0; i<A.nA; i++) 
    for (j=0; j<lda; j++) A1[i][j] = 0;

  for (i=0; i<A.nA; i++)
  {  A1[i][0]=-1; 
     A1[i][ku-1]=-1;
     A1[i][ku]=4;
     A1[i][ku+1]=-1;
     A1[i][lda-1]=-1;}   
  
  clock0 = clock();
  
  for (k=0; k<nit; k++)
     dgbmv_(&trans,&nn,&nn,&ku,&ku,&alpha,&(A1[0][0]),&lda,x,&inc,&beta,y,&inc);

  clock1 = clock();

  if (N<6) { for (i=0; i<A.nA; i++)  printf(" %e\n",y[i]);};
  printf("\n Blas CPU-time :  %e miliseconds\n",
    (1.0e+3*(double)(clock1 - clock0))/CLOCKS_PER_SEC);

  return 0;
    
}