Rosetta Code

Matrix multiplication: Difference between revisions

Page Discussion

Matrix multiplication (view source)

Revision as of 20:55, 5 May 2022

100 bytes removed ,  2 years ago
→‎{{header|C}}: Faster heapless implementation
(→‎JS ES6: Added a comment to the dotProduct function.)
(→‎{{header|C}}: Faster heapless implementation)
Line 1,085:
 
=={{header|C}}==
For performance critical work involving matrices, especially large or sparse ones, always consider using an established library such as BLAS first.
<lang c>#include <stdio.h>
#include <stdlib.h>
 
#define MAT_ELEM(rows,cols,r,c) (r*cols+c)
/* Make the data structure self-contained. Element at row i and col j
is x[i * w + j]. More often than not, though, you might want
to represent a matrix some other way */
typedef struct { int h, w; double *x;} matrix_t, *matrix;
 
//Improve performance by assuming output matrices do not overlap with
inline double dot(double *a, double *b, int len, int step)
//input matrices. If this is C++, use the __restrict extension instead
#ifdef __cplusplus
typedef double * const __restrict MAT_OUT_t;
#else
typedef double * const restrict MAT_OUT_t;
#endif
typedef const double * const MAT_IN_t;
 
static inline void mat_mul(
const int m,
const int n,
const int p,
MAT_IN_t a,
MAT_IN_t b,
MAT_OUT_t c)
{
for (int row=0; row<m; row++) {
double r = 0;
for (int col=0; col<p; col++) {
while (len--) {
c[MAT_ELEM(m,p,row,col)] = 0;
r += *a++ * *b;
for (int i = 0; i < n; i++) c[MAT_ELEM(m,p,row,col)] += a[MAT_ELEM(m,n,row,i)]*b[MAT_ELEM(n,p,i,col)];
b += step;
}
}
}
return r;
}
 
static inline void mat_show(
matrix mat_new(int h, int w)
const int m,
const int p,
MAT_IN_t c)
{
for (int row=0; row<m;row++) {
matrix r = malloc(sizeof(matrix_t) + sizeof(double) * w * h);
for (int col=0; col<p;col++) {
r->h = h, r->w = w;
printf("\t%7.3f", c[MAT_ELEM(m,p,row,col)]);
r->x = (double*)(r + 1);
}
return r;
putchar('\n');
}
}
 
int main(void)
matrix mat_mul(matrix a, matrix b)
{
double a[4*4] = {1, 1, 1, 1,
matrix r;
2, 4, 8, 16,
double *p, *pa;
3, 9, 27, 81,
int i, j;
4, 16, 64, 256};
if (a->w != b->h) return 0;
 
double b[4*3] = { 4.0, -3.0, 4.0/3,
r = mat_new(a->h, b->w);
-13.0/3, 19.0/4, -7.0/3,
p = r->x;
for (pa = a->x, i = 0; i < a->h; i++ 3.0/2, pa += a->w)2.0, 7.0/6,
-1.0/6, 1.0/4, -1.0/6};
for (j = 0; j < b->w; j++)
*p++ = dot(pa, b->x + j, a->w, b->w);
return r;
}
 
double c[4*3] = {0};
void mat_show(matrix a)
{
mat_mul(4,4,3,a,b,c);
int i, j;
mat_show(4,3,c);
double *p = a->x;
return 0;
for (i = 0; i < a->h; i++, putchar('\n'))
for (j = 0; j < a->w; j++)
printf("\t%7.3f", *p++);
putchar('\n');
}
 
</lang>
int main()
{
double da[] = { 1, 1, 1, 1,
2, 4, 8, 16,
3, 9, 27, 81,
4,16, 64, 256 };
double db[] = { 4.0, -3.0, 4.0/3,
-13.0/3, 19.0/4, -7.0/3,
3.0/2, -2.0, 7.0/6,
-1.0/6, 1.0/4, -1.0/6};
 
matrix_t a = { 4, 4, da }, b = { 4, 3, db };
matrix c = mat_mul(&a, &b);
 
/* mat_show(&a), mat_show(&b); */
mat_show(c);
/* free(c) */
return 0;
}</lang>
 
=={{header|C sharp|C#}}==
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/49275"