# Some BLAS Level 1 Functions

• For the sake of simplicity we will only consider matrix/vector elements of type double in this class.

• All vectors must be allocated before the function gets called.

• An operator description like $$x \rightarrow y$$ is read as vector $$y$$ gets overwritten with vector $$y$$.

• Furthermore:

• lower case letters indicate vectors, upper case letters matrices.

• it is assumed that different letters refer to non-overlapping vectors/matrices, i.e. they have no elements in common. Otherwise the result might be undefined.

All BLAS function names get an additional prefix that depends on the matrix/vector element type:

 Prefix Meaning Element Type in C s single precision float d double precision double c complex single precision float _Complex z complex double precision double _Complex

So functions scopy, dcopy, ccopy and zcopy provide the same functionality for vectors where elements are of type float, double, float _Complex and double _Complex respectively.

## Short Description of some BLAS Level 1 Functions

### copy: Copy Vectors

 Operation Signature $$x \rightarrow y$$ void dcopy(size_t n, const double *x, ptrdiff_t incX, double *y, ptrdiff_t incY); 

### swap: Interchanging Vectors

 Operation Signature $$x \leftrightarrow y$$ void dswap(size_t n, double *x, ptrdiff_t incX, double *y, ptrdiff_t incY); 

### axpy: Adding Vectors (Alpha x Plus y)

 Operation Signature $$\alpha x+y\rightarrow y$$ void daxpy(size_t n, double alpha, const double *x, ptrdiff_t incX, double *y, ptrdiff_t incY); 
• In the special case $$\alpha=0$$ the function returns immediately (NOP case).

### scal: Scaling Vectors

 Operation Signature $$\alpha x\rightarrow x$$ void dscal(size_t n, double alpha, double *x, ptrdiff_t incX); 
• In the special case $$\alpha=1$$ the function returns immediately (NOP case).

• In the special case $$\alpha=0$$ must not multiply elements of $$x$$ with zero. Instead $$x$$ gets overwritten with zeros. Hence, in this case we allow that $$x$$ contains NaN values.

### dot: dot product

 Operation Signature $$x^T y \rightarrow alpha$$ double ddot(size_t n, const double *x, ptrdiff_t incX, const double *y, ptrdiff_t incY); 

## Exercise

Complete the following code by

• implementing the BLAS Level 1 functions and

• doing some (very) simple test calls in function main

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>

void
initMatrix(size_t m, size_t n,
double *A,
ptrdiff_t incRow, ptrdiff_t incCol)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
A[i*incRow+j*incCol] = i*n+j+1;
}
}
}

void
printMatrix(size_t m, size_t n,
const double *A,
ptrdiff_t incRow, ptrdiff_t incCol)
{
for (size_t i=0; i<m; ++i) {
for (size_t j=0; j<n; ++j) {
printf("%10.3lf", A[i*incRow+j*incCol]);
}
printf("\n");
}
printf("\n");
}

//------------------------------------------------------------------------------

void
dcopy(size_t n,
const double *x, ptrdiff_t incX,
double *y, ptrdiff_t incY)
{
// *** Your code here ***
}

void
dswap(size_t n,
double *x, ptrdiff_t incX,
double *y, ptrdiff_t incY)
{
// *** Your code here ***
}

void
daxpy(size_t n, double alpha,
const double *x, ptrdiff_t incX,
double *y, ptrdiff_t incY)
{
// *** Your code here ***
}

void
dscal(size_t n, double alpha,
double *x, ptrdiff_t incX)
{
// *** Your code here ***
}

double
ddot(size_t n,
const double *x, ptrdiff_t incX,
const double *y, ptrdiff_t incY)
{
// *** Your code here ***
}

int
main()
{
// Allocate matrix A with dimension m x n
size_t m = 5;
size_t n = 8;

size_t incRow = n;
size_t incCol = 1;

double *A = (double *) malloc(m*n*sizeof(double));

// Initialize matrix A
initMatrix(m, n, A, incRow, incCol);

// Print initial matrix A
printf("A=\n");
printMatrix(m, n, A, incRow, incCol);

// Some simple tests for the BLAS Level 1 functions

// 1) dcopy
printf("A(:,2) -> A(:, 3)\n");
// *** Your code here ***
printMatrix(m, n, A, incRow, incCol);

// 2) dswap
printf("A(3,:) <-> A(1,:)\n");
// *** Your code here ***
printMatrix(m, n, A, incRow, incCol);

// 3) daxpy
printf("2.5*A(0,:) + A(2,:) -> A(2,:)\n");
// *** Your code here ***
printMatrix(m, n, A, incRow, incCol);

// 4) dscal
printf("0.5*A(:,5) -> A(:,5)\n");
// *** Your code here ***
printMatrix(m, n, A, incRow, incCol);

// Deallocate matrix A
free(A);
}