driver_gemv.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#define dabs( x ) ( (x) < 0 ? -(x) : x )

double FLA_Clock();      // This is a routine for extracting elapsed
			 // time borrowed from the libflame library

/* MaxAbsDiff computes the maximum absolute difference over all
   corresponding elements of two matrices */
double MaxAbsDiff( int, int, double *, int, double *, int );

/* RandomMatrix overwrites a matrix with random values */
void RandomMatrix( int, int, double *, int );

/* Prototype for BLAS matrix-matrix multiplication routine (which we will 
   use for the reference implementation */
void dgemv_( char *,                         // transA
	     int *, int *       ,            // m, n
	     double *, double *, int *,      // alpha, A, ldA
	               double *, int *,      //        x, incx
	     double *, double *, int * );    // beta,  y, incy

/* MyGemv is a common interface to all the implementations we will 
   develop so we don't have to keep rewriting this driver routine. */
void MyGemv( int, int, double *, int, double *, int, double *, int );

int main(int argc, char *argv[])
{
  int
    m, n,
    ldA, 
    size, first, last, inc,
    i_one=1,
    AllCorrect, 
    err;

  double
    d_one = 1.0,
    diff, maxdiff;

  double
    *A, *x, *y, *yold, *yref;

  /* Trials for matrix sizes m=n=first to last in increments
     of inc will be performed.  (Actually, we are going to go from
     largest to smallest since this seems to give more reliable 
     timings.  */
  printf( "%% enter first, last, inc:" );
  err = scanf( "%d%d%d", &first, &last, &inc );

  /* Adjust first and last so that they are multiples of inc */
  last = ( last / inc ) * inc;
  first = ( first / inc ) * inc;
  first = ( first == 0 ? inc : first );
  
  printf( "%% %d %d %d \n", first, last, inc );

  /* AllCorrect keeps track of whether all the trials were correct so
     far */
 
  AllCorrect = 1;

  for ( size=last; size>= first; size-=inc ){
    /* We will only test cases where all three matrices are square.
     obviously, this is not very complete. */
    m = n = size;
    ldA = size;

    /* Allocate space for the matrix and vectors. 
       A will be the address where A is stored.   
       x will be the address where x is stored.   
       y will be the address where y is stored.   

       Now, we will compute y = A x + y with via routine MyGemv
       and also with a reference implementation.  Therefore, we will
       utilize two more arrays:
 
       yold will be the address where the original vector y is
       stored.  

       yref will be the address where the result of computing y = A x
       + y computed with the reference implementation will be stored.
    */

    A = ( double * ) malloc( ldA * n * sizeof( double ) );
    x = ( double * ) malloc( n * sizeof( double ) );
    y = ( double * ) malloc( m * sizeof( double ) );
    yold = ( double * ) malloc( m * sizeof( double ) );
    yref = ( double * ) malloc( m * sizeof( double ) );

    /* Generate random matrix A */
    RandomMatrix( m, n, A, ldA );

    /* Generate random vector x */
    RandomMatrix( n, 1, x, n );

    /* Generate random vector y */
    RandomMatrix( m, 1, yold, m );
    
    /* Time reference implementation provided by the BLAS library
       routine dgemv (double precision general matrix-vector
       multiplicationn */
      
    /* Copy matrix Cold to Cref */
    memcpy( yref, yold, m * sizeof( double ) );
    
    /* Compute yref = A x + yref */
    dgemv_( "No transpose", 
	    &m, &n, 
	    &d_one, A, &ldA,
	            x, &i_one,
	    &d_one, yref, &i_one );

    /* Time MyGemv */

    /* Copy vector yold to y */
    memcpy( y, yold, m * sizeof( double ) );
    
    /* Compute y = A x + y */
    MyGemv( m, n, A, ldA, x, 1, y, 1 );

    diff = MaxAbsDiff( m, 1, y, m, yref, m );
    maxdiff = ( diff > maxdiff ? diff : maxdiff );

    if ( maxdiff > 1.0e-8 )
      /* Looks fishy... */
      AllCorrect = 0;

    /* Free the buffers */
    free( A );
    free( x );
    free( y );
    free( yold );
    free( yref );
  }

  if ( AllCorrect )
    printf("It appears all is well\n");
  else
    printf("At least one test is suspicious\n");

  exit( 0 );
}