Content |
First BLAS Level 2 Function
-
In a first approach we will implement the BLAS Level 2 functions by using Level 1 functions as building blocks.
-
Each function should distinguish the cases that a matrix is row or column major.
ger: General Matrix Rank 1 Update
Operation |
Signature |
\(\alpha xy^T+A\rightarrow A\) |
void dger(size_t m, size_t n, double alpha, const double *x, ptrdiff_t incX, const double *y, ptrdiff_t incY, double *A, ptrdiff_t incRowA, ptrdiff_t incColA); |
In the special case \(\alpha = 0\) or \(m=0\) or \(n=0\) this is a NOP.
Exercise
Complete the following code
#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"); } //-- BLAS Level 1 -------------------------------------------------------------- void dcopy(size_t n, const double *x, ptrdiff_t incX, double *y, ptrdiff_t incY) { for (size_t i=0; i<n; ++i) { y[i*incY] = x[i*incX]; } } void dswap(size_t n, double *x, ptrdiff_t incX, double *y, ptrdiff_t incY) { for (size_t i=0; i<n; ++i) { double tmp = x[i*incX]; x[i*incX] = y[i*incY]; y[i*incY] = tmp; } } void daxpy(size_t n, double alpha, const double *x, ptrdiff_t incX, double *y, ptrdiff_t incY) { if (alpha==0) { return; } for (size_t i=0; i<n; ++i) { y[i*incY] += alpha*x[i*incX]; } } void dscal(size_t n, double alpha, double *x, ptrdiff_t incX) { if (alpha==1) { return; } if (alpha==0) { for (size_t i=0; i<n; ++i) { x[i*incX] = 0; } } else { for (size_t i=0; i<n; ++i) { x[i*incX] *= alpha; } } } double ddot(size_t n, const double *x, ptrdiff_t incX, const double *y, ptrdiff_t incY) { double alpha = 0; for (size_t i=0; i<n; ++i) { alpha += x[i*incX]*y[i*incY]; } return alpha; } //-- BLAS Level 2 -------------------------------------------------------------- void dger(size_t m, size_t n, double alpha, const double *x, ptrdiff_t incX, const double *y, ptrdiff_t incY, double *A, ptrdiff_t incRowA, ptrdiff_t incColA) { // ***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); // Matrix A refers to A_(1:m-1, 1:n-1) size_t m = m_ - 1; size_t n = n_ - 1; double *A = 0; /* ***Fix me *** */ ptrdiff_t incRowA = 0; /* ***Fix me *** */ ptrdiff_t incColA = 0; /* ***Fix me *** */ // Vector x refers to A_(0, 1:n-1) double *x = 0; /* ***Fix me *** */ ptrdiff_t incX = 0; /* ***Fix me *** */ // Vector y refers to A_(1:m-1, 0) double *y = 0; /* ***Fix me *** */ ptrdiff_t incY = 0; /* ***Fix me *** */ // print A, x, y (for saving space: display vectors as row vectors) printf("A=\n"); printMatrix(m, n, A, incRowA, incColA); printf("x=\n"); printMatrix(1, n, x, 1, incX); printf("y=\n"); printMatrix(1, m, y, 1, incY); // Some simple tests for the BLAS Level 2 functions // 1) dger printf("1.5*y * x^T + A -> A\n"); // ***Your code here*** printMatrix(m, n, A, incRow, incCol); // 2) dger printf("1.5*x * y^T + A^T -> A^T\n"); // ***Your code here*** printMatrix(m, n, A, incRow, incCol); // Deallocate matrix A_ free(A_); }