Call Today 716.688.4675

Basic Linear Algebra Subprograms

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:

Dot Product Definition

To implement this in C, we use the cdotc function:

CDOTC function

So for example, let us defi ne:

DotEx1

The dot product is then given by:

DotEx2

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:

Matrix Vector Definition

To perform this operation using CBLAS, we use the cgemv function:

CGEMVEX

So, for example, let us defi ne:

Matrix Vector Ex1

So then we have:

Matrix Vector Ex2

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:

MatMatProdDef

To perform this operation using CBLAS, we use the cgemm function:

CGEMMDEF

So, for example, defi ne:

CGEMMEX

Then we have:

CGEMMEX1

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:

CGESV

So, for example, let us de fine:

CGESV EX

Then we have:

CGESV EX1

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[0] = 1.0 + 1.0*I;
a[1] = 2.0 – 1.0*I;

b[0] = 3.0 – 2.0*I;
b[1] = 1.0 + 1.0*I;

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

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===//
printf(“ACML Answer b:\n”);
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===//
printf(“CBLAS Answer b:\n”);
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);
}

VOCAL Technologies, Ltd.
520 Lee Entrance, Suite 202
Amherst New York 14228
Phone: +1-716-688-4675
Fax: +1-716-639-0713
Email: sales@vocal.com