Basic Linear Algebra Subprograms (BLAS) are a library of low-level linear algebra subroutines including dot products, matrix vector multiplication, and matrix matrix multiplication. Examples are provided for the C language extension CBLAS using single precision complex arithmetic to show how to perform dot products, matrix vector multiplication, matrix matrix multiplication, as well as solve a system of linear equations. For ease of implementation, complex.h is used where possible. In addition, examples using the optimized ACML (AMD Core Math Library) are provided.

## Dot Product

The dot product of two complex vectors is defi ned as: To implement this in C, we use the cdotc function: So for example, let us defi ne: The dot product is then given by: ## Matrix Vector Product

Suppose we have a matrix A of size NxM, in other words, it has N rows and M columns. The vector x that it multiplies, must be of size M, in other words, it must have M elements. The resulting vector b will then be of size N. The calculation procedes as follows: To perform this operation using CBLAS, we use the cgemv function: So, for example, let us defi ne: So then we have: ## Matrix Matrix Product

Suppose we have a matrix A of size NxM, and another matrix B of size MxP. The result of their multiplication, C, is then of size NxP. The calculation proceeds as follows: To perform this operation using CBLAS, we use the cgemm function: So, for example, defi ne: Then we have: ## Solving Systems of Linear Equations

Suppose we have a matrix A of size NxN, another matrix B of size NxNRHS, and we wish to compute the solution of AX = B for the matrix X of size NxNRHS. We can so do by fi rst performing a LU decomposition with partial pivoting. Thus, we express A as A = P · L · U where P is the permutation matrix, L is the lower triangular matrix, and U is the upper triangular matrix.

Solving a system of linear equations in this manner cannot actually be done with the basic CBLAS routines. However, the ACML does provide a means to perform this operation, namely, the XGESV functions, similar to the LAPACK extension to CBLAS. To nd the solution, we use the function cgesv: So, for example, let us de fine: Then we have: ## Example 1: CDOTC

#include
#include
#include
#include “acml.h”
#include “cblas.h”

#define N 2

int main(int argc, char *argv[]){

int inca = 1, incb = 1;
complex c;
float _Complex a[N], b[N], r;
void *aptr, *bptr, *rptr;
aptr = &a, bptr = &b; rptr = &r;

printf(“\n\\======== An example to illustrate the dot product of vectors (a|b)=c =======\\\n\n”);

a = 1.0 + 1.0*I;
a = 2.0 – 1.0*I;

b = 3.0 – 2.0*I;
b = 1.0 + 1.0*I;

printf(“a = ( %6.2f, %6.2f) ( %6.2f, %6.2f) \n”,creal(a),cimag(a),creal(a),cimag(a));
printf(“b = ( %6.2f, %6.2f) ( %6.2f, %6.2f) \n”,creal(b),cimag(b),creal(b),cimag(b));

printf(“\n”);
cblas_cdotc_sub(N, aptr, 1, bptr, 1, rptr);
printf(“The cblas dot product (a|b) is: ( %6.2f, %6.2f) \n”, creal(r), cimag(r) );
c = cdotc(N, a, inca, b, incb);
printf(“The acml dot product (a|b) is: ( %6.2f, %6.2f) \n”, c.real, c.imag );
printf(“\n”);

return(0);
}

## Example 2: CGEMV

#include
#include
#include
#include
#include “acml.h”
#include “cblas.h”

#define NUMROWS 3
#define NUMCOLS 2

#define A(I,J) a[I + J*NUMROWS]

int main(int argc, char *argv[]){

int i,j;
float _Complex a[NUMROWS*NUMCOLS], x[NUMCOLS], b[NUMROWS];
float _Complex alpha = 1.0 + 0.0*I, beta = 0.0 + 0.0*I;
void *AComplex, *alphaComplex, *XComplex, *betaComplex, *BComplex;
AComplex = &a; alphaComplex = α XComplex = &x; betaComplex = β BComplex = &b;

printf(“\n\\======= An example to illustrate matrix vector multiplication Ax=b =======\\\n\n”);

//===Fill the Data===//
for (i=0; i<NUMROWS*NUMCOLS; i++) a[i] = i + i*I;
for (i=0; i<NUMCOLS; i++) x[i] = (NUMCOLS-i) + (NUMCOLS-i)*I;
memset(b,0,NUMROWS*sizeof(float _Complex));

//===Perform thhe ACML Matrix Vector Multiplication===//
cgemv(‘N’, NUMROWS, NUMCOLS, &alpha, a, NUMROWS, x, 1, &beta, b, 1);

//===Print the Data===//
printf(“Matrix A:\n”);
for (i = 0; i < NUMROWS; i++){
for (j = 0; j < NUMCOLS; j++){
printf(“(%3.2f,%3.2f)”, creal(A(i,j)), cimag(A(i,j)));
}
printf(“\n”);
}
printf(“\n”);

printf(“Vector x:\n”);
for (i = 0; i < NUMCOLS; i++){
printf(“(%3.2f,%3.2f) \n”, creal(x[i]), cimag(x[i]));
}
printf(“\n”);

//===Print the solution===//
for (i = 0; i < NUMROWS; i++){
printf(“(%3.2f,%3.2f) \n”, creal(b[i]), cimag(b[i]));
}
printf(“\n”);

//===Reinitialize the solution===//
memset(b,0,NUMROWS*sizeof(float _Complex));

//===Perform CBLAS Matrix Vector Multiplication===//
cblas_cgemv(CblasColMajor,CblasNoTrans, NUMROWS, NUMCOLS,
alphaComplex, AComplex, NUMROWS, XComplex, 1,
betaComplex, BComplex, 1);

//===Print the CBLAS solution===//
for (i = 0; i < NUMROWS; i++){
printf(“(%3.2f,%3.2f) \n”, creal(b[i]), cimag(b[i]));
}

return(0);
}

## Example 3: CGEMM

#include
#include
#include
#include
#include “acml.h”
#include “cblas.h”

#define NUMROWS_A 3
#define NUMCOLS_A 2
#define NUMROWS_B 2
#define NUMCOLS_B 3

#define A(I,J) a[I +J*NUMROWS_A]
#define B(I,J) b[I +J*NUMROWS_B]
#define C(I,J) c[I +J*NUMROWS_A]

int main(int argc, char *argv[]){
int i,j,count;
float _Complex a[NUMROWS_A*NUMCOLS_A], b[NUMROWS_B*NUMCOLS_B], c[NUMROWS_A*NUMCOLS_B];
float _Complex alpha = 1.0 + 0.0*I, beta = 0.0 + 0.0*I;
void *AComplex, *alphaComplex, *BComplex, *betaComplex, *CComplex;
AComplex = &a; alphaComplex = α BComplex = &b; betaComplex = β CComplex = &c;

printf(“\n\\==== An example to illustrate matrix matrix multiplication AB = C =====\\\n\n”);

//===Initialize the Data===//
for (i=0; i<NUMROWS_A*NUMCOLS_A; i++) a[i] = i + i*I;
count = NUMROWS_B*NUMCOLS_B-1;
for (i=0; i<NUMROWS_B*NUMCOLS_B; i++) b[i] = (1.0 + 1.0*I)*(count–);
memset(c,0,NUMROWS_A*NUMCOLS_B*sizeof(float _Complex));

//===Perform the CBLAS Matrix Matrix Multiplication===//
cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
NUMROWS_A, NUMCOLS_B, NUMCOLS_A, alphaComplex,
AComplex, NUMROWS_A, BComplex, NUMROWS_B, betaComplex,
CComplex, NUMROWS_A);

//===Print the Data===//
printf(“Matrix A:\n”);
for (i = 0; i < NUMROWS_A; i++){
for (j = 0; j < NUMCOLS_A; j++){
printf(“(%3.2f, %3.2f) “, creal(A(i,j)), cimag(A(i,j)) );
}
printf(“\n”);
}
printf(“\n”);
printf(“Matrix B:\n”);
for (i = 0; i < NUMROWS_B; i++){
for (j = 0; j < NUMCOLS_B; j++){
printf(“(%3.2f, %3.2f) “, creal(B(i,j)), cimag(B(i,j)) );
}
printf(“\n”);
}
//===Print the CBLAS Result===//
printf(“\nCBLAS Result:\n”);
for (i = 0; i < NUMROWS_A; i++){
for (j = 0; j < NUMCOLS_B; j++){
printf(“(%3.2f, %3.2f) “, creal(C(i,j)), cimag(C(i,j)) );
}
printf(“\n”);
}

//===Reinitialize the Solution===//
memset(c,0,NUMROWS_A*NUMCOLS_B*sizeof(float _Complex));

//===Perform the ACML Matrix Matrix Multiplication===//
cgemm(‘N’, ‘N’, NUMROWS_A, NUMCOLS_B, NUMCOLS_A, &alpha,
a, NUMROWS_A, b, NUMROWS_B, &beta, c, NUMROWS_A);

//===Print the ACML Result===//
printf(“\nACML Result:\n”);
for (i = 0; i < NUMROWS_A; i++){
for (j = 0; j < NUMCOLS_B; j++){
printf(“(%3.2f, %3.2f) “, creal(C(i,j)), cimag(C(i,j)) );
}
printf(“\n”);
}
return(0);
}

## Example 4: CGESV

#include
#include
#include
#include “acml.h”

#define NUMCOLS 3
#define NUMROWS 3
#define NUMRTHS 1

#define A(I,J) a[I + NUMROWS*J]
#define B(I,J) b[I + NUMROWS*J]

int main(int argc, char *argv[]){

printf(“\n\\====== An example illustrating how to compute the solution of Ax=b ======\\ \n\n”);

int i,j,info;
int pivot[NUMROWS];
float _Complex a[NUMROWS*NUMCOLS], b[NUMROWS*NUMRTHS];

//===Fill the Matrix A===//
A(0,0) = 2.0 + 1.0*I; A(0,1) = 0.0 + 0.0*I; A(0,2) = 1.0 + 0.0*I;
A(1,0) = 3.0 + 0.0*I; A(1,1) = 0.0 + 1.0*I; A(1,2) = 0.0 + 0.0*I;
A(2,0) = 5.0 + 0.0*I; A(2,1) = 1.0 + 0.0*I; A(2,2) = 1.0 + 1.0*I;

//===Fill the Right Hand Side b==//
B(0,0) = 3.0 + 3.0*I;
B(1,0) = 2.0 + 2.0*I;
B(2,0) = 1.0 + 1.0*I;

//===Print the Data before Solving===//
printf(“Matrix A:\n”);
for (i = 0; i < 3; i++){
for (j = 0; j < 3; j++){
printf(” (%3.2f,%3.2f)”, creal(A(i,j)), cimag(A(i,j)));
}
printf(“\n”);
}
printf(“\n”);
printf(“Right Hand Side b:\n”);
for (i = 0; i < 3; i++){
for (j = 0; j < 1; j++){
printf(” (%3.2f,%3.2f)”, creal(B(i,j)), cimag(B(i,j)));
}
printf(“\n”);
}
printf(“\n”);

//===Compute solution of Ax = b ===//
cgesv(NUMROWS, NUMRTHS, a, NUMROWS, pivot, b, NUMROWS, &info);

//===Print Solution===//
printf(“Solution x:\n”);
for (i = 0; i < 3; i++){
for (j = 0; j < 1; j++){
printf(” (%3.2f,%3.2f)”, creal(B(i,j)), cimag(B(i,j)));
}
printf(“\n”);
}

return(0);
}