Matrix multiplication: Difference between revisions

Content added Content deleted

Inline

@@ Line 342: / Line 342: @@
 =={{header|C}}==
+For performance critical work involving matrices, especially large or sparse ones, always consider using an established library such as BLAS first.
-{{works with|gcc|<nowiki>4.1.2 20070925 (Red Hat 4.1.2-27) Options: gcc -std=gnu99</nowiki>}}
 <lang c>#include <stdio.h>
+#include <stdlib.h>
-#define dim 4 /* fixed length square matrices */
-const int SLICE=0; /* coder hints */
-typedef double field_t; /* field_t type is long float */
-typedef field_t vec_t[dim];
-typedef field_t *ref_vec_t; /* address of first element */
-typedef vec_t matrix_t[dim];
-typedef field_t *ref_matrix_t; /* address of first element */
-typedef const char* format;
+/* Make the data structure self-contained.  Element at row i and col j
-/* define the vector/matrix_t operators */
+   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;
-field_t v_times_v (vec_t a, vec_t b, int step_b){
+inline double dot(double *a, double *b, int len, int step)
+{
-    /* basically the dot product if step_b==1*/
+	double r = 0;
-    field_t result=0;
+	while (len--) {
-    for( int i=0; i<sizeof a; i++ )
-      result+= a[i]*b[i*step_b];
+		r += *a++ * *b;
+		b += step;
-    return result;
-  }
+	}
+	return r;
+}
+matrix mat_new(int h, int w)
-ref_vec_t v_eq_v_times_m(vec_t result, vec_t a, matrix_t b){
+{
-    for( int j=0; j<sizeof b; j++ )
-      result[j]=v_times_v(a,&b[SLICE][j],sizeof b[SLICE] / sizeof (field_t));
+	matrix r = malloc(sizeof(matrix) + sizeof(double) * w * h);
+	r->h = h, r->w = w;
-    return &result[SLICE];
+	r->x = (double*)(r + 1);
-  }
+	return r;
+}
+matrix mat_mul(matrix a, matrix b)
-ref_matrix_t m_eq_m_times_m (matrix_t result, matrix_t a, matrix_t b){
+{
-    for( int k=0; k<sizeof result; k++ )
+	matrix r;
-      v_eq_v_times_m(&result[k][SLICE],&a[k][SLICE],b);
+	double *p, *pa;
-    return &result[SLICE][SLICE];
+	int i, j;
-  }
+	if (a->w != b->h) return 0;
+	r = mat_new(a->h, b->w);
-/* Some sample matrices to test */
+	p = r->x;
-matrix_t a={{1,  1,  1,   1}, /* matrix_t A */
-            {2,  4,  8,  16},
+	for (pa = a->x, i = 0; i < a->h; i++, pa += a->w)
+		for (j = 0; j < b->w; j++)
-            {3,  9, 27,  81},
+			*p++ = dot(pa, b->x + j, a->w, b->w);
-            {4, 16, 64, 256}};
+	return r;
+}
+void mat_show(matrix a)
-matrix_t b={{  4.0  , -3.0  ,  4.0/3,  -1.0/4 }, /* matrix_t B */
+{
-            {-13.0/3, 19.0/4, -7.0/3,  11.0/24},
+	int i, j;
-            {  3.0/2, -2.0  ,  7.0/6,  -1.0/4 },
+	double *p = a->x;
-            { -1.0/6,  1.0/4, -1.0/6,   1.0/24}};
+	for (i = 0; i < a->h; i++, putchar('\n'))
+		for (j = 0; j < a->w; j++)
+			printf("\t%7.3f", *p++);
+	putchar('\n');
+}
-int main(){
+int main()
+{
-  matrix_t prod;
+	double da[] = {	1, 1,  1,   1,
-  m_eq_m_times_m(prod,a,b); /* actual multiplication example of A x B */
+, 4,  8,  16,
+, 9, 27,  81,
+,16, 64, 256	};
+	double db[] = {     4.0,   -3.0,  4.0/3,
+			-13.0/3, 19.0/4, -7.0/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 };
-  #define field_fmt "%6.2f" /* width of 6, with no '+' sign, 2 decimals */
+	matrix c = mat_mul(&a, &b);
-  #define vec_fmt "{"field_fmt","field_fmt","field_fmt","field_fmt"}"
-  #define matrix_fmt " {"vec_fmt",\n  "vec_fmt",\n  "vec_fmt",\n  "vec_fmt"};"
-  format result_fmt = " Product of a and b: \n"matrix_fmt"\n";
+	/* mat_show(&a), mat_show(&b); */
-  /* finally print the result */
+	mat_show(c);
-  vprintf(result_fmt,(void*)&prod);
+	/* free(c) */
+	return 0;
 }</lang>
-Output:
- Product of a and b:
- {{  1.00,  0.00, -0.00, -0.00},
-  {  0.00,  1.00, -0.00, -0.00},
-  {  0.00,  0.00,  1.00, -0.00},
-  {  0.00,  0.00,  0.00,  1.00}};
 =={{header|C sharp|C#}}==