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[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);
}