I wrote this simple C file years ago, then when I needed to remember how to multiply vectors and matrices using Fortran routines in CBLAS/ LAPACK / ATLAS I went looking for it again.
A simple reference is here http://www.seehuhn.de/pages/linear
The biggest mistake with DGEMV is that incrementX and incrementY is the column of the vector, or 1 and NOT the rank of the vector unless your vector is many parallel vectors.
Many people put m or n, the rank of the matrix, into the vector values incrementX and incrementY but you should use 1 for both unless you have a row of vectors in an array.
BLAS was made to do multiple dimension linear algebra computations very efficiently. So it's laid out counter-intuitively for people just trying to multiply two small examples.
Compile this simple program to see how the linear algebra variables are handled for yourself.
A simple reference is here http://www.seehuhn.de/pages/linear
The biggest mistake with DGEMV is that incrementX and incrementY is the column of the vector, or 1 and NOT the rank of the vector unless your vector is many parallel vectors.
Many people put m or n, the rank of the matrix, into the vector values incrementX and incrementY but you should use 1 for both unless you have a row of vectors in an array.
BLAS was made to do multiple dimension linear algebra computations very efficiently. So it's laid out counter-intuitively for people just trying to multiply two small examples.
Compile this simple program to see how the linear algebra variables are handled for yourself.
#include <stdio.h>
#include <cblas.h>
#include <clapack.h>
#define DE
int main (void)
{
// testing cblas function calls
// this test file shows how to use the CBLAS and LAPACK function calls
// The solutions here were confirmed with MATLAB
// matrices
// R1C1 R1C2 R1C3
double a[] = { 0.11, 0.12, 0.13,
0.21, 0.22, 0.23 };
double b[] = { 1011, 1012,
1021, 1022,
1031, 1032 };
double c[] = { 0.00, 0.00,
0.00, 0.00 };
double d[] = { 2.00, 0.00,
0.00, 2.00 };
double e[] = { 0.10, 0.20,
0.30, 0.40 };
double f[] = { 1.0, 0.0,
0.0, 1.0 };
double g[] = { 1.00, 3.00,
2.00, 4.00 };
double i[] = { 1.00, 0.00,
0.00, 1.00 };
double j[] = { 2.00, 3.00,
3.00, 2.00 };
// vectors for experimentation
double y[] = { 1.0, 1.0};
double z[] = { 1.0, 1.0};
// pivot table for 2 by 2
int p[]= { 1.0, 1.0 };
// The BLAS 2D matrices access all values along the 1st column in sequence from the array
// and then the next column etc. until the (n-1)th column of the matrix. It uses the
// indices values, I believe, to stop reading in matrix values.
// Columns x Rows so X(i,j) means column i and
// matrix scaling
double alpha = 1.0;
double beta = 0.0;
printf("testing cblas function calls \n");
#ifdef DE
printf ("[ %g, %g, %g\n", a[0], a[1], a[2]);
printf (" %g, %g, %g ]\n*\n", a[3], a[4], a[5]);
printf ("[ %g, %g\n", b[0], b[1]);
printf (" %g, %g \n", b[2], b[3]);
printf (" %g, %g ]\n=\n", b[4], b[5]);
#endif
/*
cblas_dgemm: Compute C = alpha*AB + beta*C
NOTES FOR DGEMM:
Matrix A is returned altered from the dgemm function
so it should be sent as a copy if you wish to keep the original.
A is returned with the proper LU factorization for the original matrix
In this example, A [2x3] * B [3x2] = C [2x2]
*/
printf (" matrix A before multiplication \n");
printf ("[ %g, %g, %g\n", a[0], a[1], a[2]);
printf (" %g, %g, %g ]\n*\n", a[3], a[4], a[5]);
cblas_dgemm ( /* which leading dimension to consider row or col - this case row */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasNoTrans,
/* m , n = matrix C dimensions [m x n] */2, 2,
/* k inner dimension of A*B */3,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&a,
/* cols of A */3,
/* matrix B */(const double *)&b,
/* cols of B */2,
/* scaling factor for C */(const double)beta,
/* Matrix C */(double *)&c,
/* rows of C */2);
printf (" A * B = \n");
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
printf (" matrix A after multiplication \n");
printf ("[ %g, %g, %g\n", a[0], a[1], a[2]);
printf (" %g, %g, %g ]\n*\n", a[3], a[4], a[5]);
/* Compute C = alpha*AB + beta*C */
cblas_dgemm ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasNoTrans,
/* m x n */2, 2,
/* k */2,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&i,
/* leading dimension of A */2,
/* matrix B */(const double *)&j,
/* leading dimension of B */2,
/* scaling factor for C */(const double)beta,
/* Matrix C */(double *)&c,
/* leading dimension of C */2);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
/* Compute C = alpha*AB + beta*C */
cblas_dgemm ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasNoTrans,
/* m x n */2, 2,
/* k */2,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&d,
/* leading dimension of A */2,
/* matrix B */(const double *)&e,
/* leading dimension of B */2,
/* scaling factor for C */(const double)beta,
/* Matrix C */(double *)&c,
/* leading dimension of C */2);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
/*
/*
Compute C = alpha*A(transpose(B)) + beta*C
This example shows how to apply a transpose to a matrix
*/
cblas_dgemm ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasTrans,
/* m x n */2, 2,
/* k */2,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&d,
/* leading dimension of A */2,
/* matrix B */(const double *)&e,
/* leading dimension of B */2,
/* scaling factor for C */(const double)beta,
/* Matrix C */(double *)&c,
/* leading dimension of C */2);
printf (" matrix multiplication with transpose B \n");
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n", c[2], c[3]);
/* Compute C = inv(D) */
/* int clapack_dgesv(const enum CBLAS_ORDER Order, const int N, const int NRHS,
double *A, const int lda, int *ipiv,
double *B, const int ldb);
A matrix is returned altered from this function
*/
printf (" matrix D before inversion \n");
printf ("[ %g, %g\n", d[0], d[1]);
printf (" %g, %g ]\n", d[2], d[3]);
// Operation required to be CblasColMajor to return the inverted matrix in the
// proper format
// Values in pivot table don't matter, they are replaced by the fcn.
clapack_dgesv ( /* which leading dimension to consider row or col */CblasColMajor,
/* size of matrices */ 2,
/* NRHS must be as great as M or n*/ 2,
/* matrix A */(double *)&d,
/* leading dimension of A */2,
/* pivot table for columns */(int *)&p,
/* matrix B */(double *)&f,
/* leading dimension of B */2);
printf (" matrix D after inversion \n");
printf ("[ %g, %g\n", d[0], d[1]);
printf (" %g, %g ]\n", d[2], d[3]);
printf (" inverted matrix \n");
printf ("[ %g, %g\n", f[0], f[1]);
printf (" %g, %g ]\n", f[2], f[3]);
// Problem with the inverse operation, it is not working for A but for transpose A.
printf ("computing A*inv(A) = \n");
/* Compute C = alpha*D(inverse(D)) + beta*C */
cblas_dgemm ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasNoTrans,
/* m x n */2, 2,
/* k */2,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&g,
/* leading dimension of A */2,
/* matrix B */(const double *)&f,
/* leading dimension of B */2,
/* scaling factor for C */(const double)beta,
/* Matrix C */(double *)&c,
/* leading dimension of C */2);
printf ("[ %g, %g\n", g[0], g[1]);
printf (" %g, %g ]\n*\n", g[2], g[3]);
printf ("[ %g, %g\n", f[0], f[1]);
printf (" %g, %g ]\n=\n", f[2], f[3]);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n\n", c[2], c[3]);
printf ("computing C= betaC \n");
/* Compute C = alpha*D(inverse(D)) + beta*C */
cblas_dgemm ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* Transpose B or not */ CblasNoTrans,
/* m x n */2, 2,
/* k */2,
/* scaling factor for A */(const double)beta,
/* matrix A */(const double *)&a,
/* leading dimension of A */2,
/* matrix B */(const double *)&b,
/* leading dimension of B */2,
/* scaling factor for C */(const double)alpha*1.5,
/* Matrix C */(double *)&c,
/* leading dimension of C */2);
printf ("[ %g, %g\n", c[0], c[1]);
printf (" %g, %g ]\n\n", c[2], c[3]);
i[0] = 0.5;
i[1] = 0.0;
i[2] = 0.0;
i[3] = 0.5;
y[0] = 1.3;
y[1] = 2.5;
z[0] = y[0];
z[1] = y[1];
printf ("x = [ %g, %g ]\n", z[0], z[1]);
printf ("A =\n");
printf ("[ %g, %g\n", i[0], i[1]);
printf (" %g, %g ]\n\n", i[2], i[3]);
printf ("computing y= alpha(Ax) + beta(y) \n");
/* Compute y = alpha*(Ax) + beta*(y) */
cblas_dgemv ( /* which leading dimension to consider row or col */CblasRowMajor,
/* Transpose A or not */ CblasNoTrans,
/* m x n */2, 2,
/* scaling factor for A */(const double)alpha,
/* matrix A */(const double *)&i,
/* leading dimension of A */2,
/* vector x */(const double *)&z,
/* increment of x */2,
/* scaling factor for y */(const double)beta,
/* vector y */(double *)&y,
/* increment of y */2);
printf ("alpha = [ %g ]\n", alpha );
printf ("beta = [ %g ]\n", beta );
printf ("x = [ %g, %g ]\n", z[0], z[1]);
printf ("A = [ %g, %g\n", i[0], i[1]);
printf (" %g, %g ]\n->\n", i[2], i[3]);
printf ("y = [ %g, %g ]\n", y[0], y[1]);
printf ("x = [ %g, %g ]\n", z[0], z[1]);
printf ("A = [ %g, %g\n", i[0], i[1]);
printf (" %g, %g ]\n->\n", i[2], i[3]);
printf ("y = [ %g, %g ]\n", y[0], y[1]);
return 0;
}
No comments:
Post a Comment