#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
//#include "mmat.h"
#include <stdbool.h>
// Matrix.h
typedef struct
{
size_t y_dim, x_dim;
float **cells;
} Matrix;
Matrix *matrix_construct(size_t y_dim, size_t x_dim); // constructs zeroed
void matrix_set_values(Matrix *dst, float *row_major_floats);
void matrix_free(Matrix *m); // frees m as well as cells
void matrix_assign(Matrix *dst, Matrix const *src); // may re-size dst
bool matrix_is_equal(Matrix const *m1, Matrix const *m2);
typedef Matrix *MatrixBinaryOperator(Matrix const *m1, Matrix const *m2);
MatrixBinaryOperator matrix_multiply;
bool matrix_operate_self(MatrixBinaryOperator *op, Matrix *m1, Matrix const *m2);
//----- end of mmat.h
Matrix *matrix_construct(size_t y_dim, size_t x_dim)
{
Matrix
*m
= malloc( sizeof *m
);
m->x_dim = x_dim;
m->y_dim = y_dim;
// assume IEEE754 floats (all bits zero = 0)
m
->cells
= malloc( y_dim
* sizeof *m
->cells
); for ( size_t ii = 0; ii < y_dim; ++ii )
m
->cells
[ii
] = calloc( x_dim
, sizeof *m
->cells
[ii
] );
return m;
}
void matrix_set_values(Matrix *dst, float *row_major_floats)
{
for (size_t y = 0; y < dst->y_dim; ++y)
for (size_t x = 0; x < dst->x_dim; ++x)
dst->cells[y][x] = *row_major_floats++;
}
bool matrix_is_equal(Matrix const *m1, Matrix const *m2)
{
if ( m1->y_dim != m2->y_dim ) return false;
if ( m1->x_dim != m2->x_dim ) return false;
for (size_t y = 0; y < m1->y_dim; ++y)
if ( memcmp(m1
->cells
[y
], m2
->cells
[y
], m1
->x_dim
* sizeof m1
->cells
[0][0]) ) return false;
return true;
}
bool matrix_operate_self(MatrixBinaryOperator *op, Matrix *m1, Matrix const *m2)
{
Matrix *result = op(m1, m2);
if ( !result ) return true;
matrix_assign(m1, result);
matrix_free(result);
return false;
}
static void matrix_deallocate(Matrix *m)
{
for ( size_t y = 0; y < m->y_dim; ++y )
m->x_dim = 0;
m->y_dim = 0;
m->cells = 0;
}
void matrix_free(Matrix *m)
{
matrix_deallocate(m);
}
void matrix_assign(Matrix *dst, Matrix const *src)
{
// note: could optimize this to not allocate if dims the same
matrix_deallocate(dst);
Matrix *new = matrix_construct(src->y_dim, src->x_dim);
for (size_t y = 0; y < src->y_dim; ++y)
memcpy(new
->cells
[y
], src
->cells
[y
], src
->x_dim
* sizeof new
->cells
[0][0]);
if ( !new
) exit(EXIT_FAILURE
); *dst = *new;
}
Matrix *matrix_multiply(Matrix const *m1, Matrix const *m2)
{
if ( m1->x_dim != m2->y_dim )
return NULL;
Matrix *new = matrix_construct(m1->y_dim, m2->x_dim);
for (size_t col = 0; col < m2->x_dim; ++col) // Each column of m2
for (size_t row = 0; row < m1->y_dim; ++row) // Each row of m1
{
// do dot-product of those two
float sum = 0;
for (size_t ii = 0; ii < m2->y_dim; ++ii)
sum += m2->cells[ii][col] * m1->cells[row][ii];
new->cells[row][col] = sum;
}
return new;
}
static void matrix_printf(char const *prompt, Matrix const *mat)
{
for (size_t y = 0; y < mat->y_dim; ++y)
{
if (y
> 0 ) printf("%20.20s ", ""); for (size_t x = 0; x < mat->x_dim; ++x)
printf("%3.3f ", mat
->cells
[y
][x
]); }
}
int main()
{
float mat1_array[3][2] = {{0, 1}, {3, 4}, {6, 7}};
float mat2_array[2][3] = {{5, 1, 2}, {3, 4, 5}};
Matrix *mat1 = matrix_construct(3, 2);
Matrix *mat2 = matrix_construct(2, 3);
matrix_set_values(mat1, (float *)&mat1_array);
matrix_set_values(mat2, (float *)&mat2_array);
Matrix *mat3 = matrix_multiply(mat1, mat2);
matrix_printf("mat1", mat1);
matrix_printf("mat2", mat2);
matrix_printf("mat3", mat3);
matrix_operate_self(matrix_multiply, mat1, mat2);
assert( matrix_is_equal
(mat1
, mat3
) );
matrix_free(mat1);
matrix_free(mat2);
matrix_free(mat3);
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPGFzc2VydC5oPgojaW5jbHVkZSA8c3RyaW5nLmg+CgovLyNpbmNsdWRlICJtbWF0LmgiCiNpbmNsdWRlIDxzdGRib29sLmg+CgovLyBNYXRyaXguaAogICAgdHlwZWRlZiBzdHJ1Y3QKICAgIHsKICAgICAgICBzaXplX3QgeV9kaW0sIHhfZGltOwogICAgICAgIGZsb2F0ICoqY2VsbHM7CiAgICB9IE1hdHJpeDsKCiAgICBNYXRyaXggKm1hdHJpeF9jb25zdHJ1Y3Qoc2l6ZV90IHlfZGltLCBzaXplX3QgeF9kaW0pOwkvLyBjb25zdHJ1Y3RzIHplcm9lZAogICAgdm9pZCBtYXRyaXhfc2V0X3ZhbHVlcyhNYXRyaXggKmRzdCwgZmxvYXQgKnJvd19tYWpvcl9mbG9hdHMpOwogICAgdm9pZCBtYXRyaXhfZnJlZShNYXRyaXggKm0pOyAgICAgLy8gZnJlZXMgbSBhcyB3ZWxsIGFzIGNlbGxzCgogICAgdm9pZCBtYXRyaXhfYXNzaWduKE1hdHJpeCAqZHN0LCBNYXRyaXggY29uc3QgKnNyYyk7ICAgLy8gbWF5IHJlLXNpemUgZHN0CiAgICBib29sIG1hdHJpeF9pc19lcXVhbChNYXRyaXggY29uc3QgKm0xLCBNYXRyaXggY29uc3QgKm0yKTsgICAKCiAgICB0eXBlZGVmIE1hdHJpeCAqTWF0cml4QmluYXJ5T3BlcmF0b3IoTWF0cml4IGNvbnN0ICptMSwgTWF0cml4IGNvbnN0ICptMik7CiAgICBNYXRyaXhCaW5hcnlPcGVyYXRvciBtYXRyaXhfbXVsdGlwbHk7CiAgICBib29sIG1hdHJpeF9vcGVyYXRlX3NlbGYoTWF0cml4QmluYXJ5T3BlcmF0b3IgKm9wLCBNYXRyaXggKm0xLCBNYXRyaXggY29uc3QgKm0yKTsKCi8vLS0tLS0gZW5kIG9mIG1tYXQuaAoKCk1hdHJpeCAqbWF0cml4X2NvbnN0cnVjdChzaXplX3QgeV9kaW0sIHNpemVfdCB4X2RpbSkKewoJTWF0cml4ICptID0gbWFsbG9jKCBzaXplb2YgKm0gKTsKCgltLT54X2RpbSA9IHhfZGltOwoJbS0+eV9kaW0gPSB5X2RpbTsKCi8vIGFzc3VtZSBJRUVFNzU0IGZsb2F0cyAoYWxsIGJpdHMgemVybyA9IDApCgltLT5jZWxscyA9IG1hbGxvYyggeV9kaW0gKiBzaXplb2YgKm0tPmNlbGxzICk7Cglmb3IgKCBzaXplX3QgaWkgPSAwOyBpaSA8IHlfZGltOyArK2lpICkKCQltLT5jZWxsc1tpaV0gPSBjYWxsb2MoIHhfZGltLCBzaXplb2YgKm0tPmNlbGxzW2lpXSApOwoKCXJldHVybiBtOwp9Cgp2b2lkIG1hdHJpeF9zZXRfdmFsdWVzKE1hdHJpeCAqZHN0LCBmbG9hdCAqcm93X21ham9yX2Zsb2F0cykKewoJZm9yIChzaXplX3QgeSA9IDA7IHkgPCBkc3QtPnlfZGltOyArK3kpCgkJZm9yIChzaXplX3QgeCA9IDA7IHggPCBkc3QtPnhfZGltOyArK3gpCgkJCWRzdC0+Y2VsbHNbeV1beF0gPSAqcm93X21ham9yX2Zsb2F0cysrOwp9Cgpib29sIG1hdHJpeF9pc19lcXVhbChNYXRyaXggY29uc3QgKm0xLCBNYXRyaXggY29uc3QgKm0yKQp7CglpZiAoIG0xLT55X2RpbSAhPSBtMi0+eV9kaW0gKSByZXR1cm4gZmFsc2U7CglpZiAoIG0xLT54X2RpbSAhPSBtMi0+eF9kaW0gKSByZXR1cm4gZmFsc2U7CgoJZm9yIChzaXplX3QgeSA9IDA7IHkgPCBtMS0+eV9kaW07ICsreSkKCQlpZiAoIG1lbWNtcChtMS0+Y2VsbHNbeV0sIG0yLT5jZWxsc1t5XSwgbTEtPnhfZGltICogc2l6ZW9mIG0xLT5jZWxsc1swXVswXSkgKQoJCQkJcmV0dXJuIGZhbHNlOwoKCXJldHVybiB0cnVlOwp9Cgpib29sIG1hdHJpeF9vcGVyYXRlX3NlbGYoTWF0cml4QmluYXJ5T3BlcmF0b3IgKm9wLCBNYXRyaXggKm0xLCBNYXRyaXggY29uc3QgKm0yKQp7CglNYXRyaXggKnJlc3VsdCA9IG9wKG0xLCBtMik7CglpZiAoICFyZXN1bHQgKSByZXR1cm4gdHJ1ZTsKCW1hdHJpeF9hc3NpZ24obTEsIHJlc3VsdCk7CgltYXRyaXhfZnJlZShyZXN1bHQpOwoJcmV0dXJuIGZhbHNlOwp9CgpzdGF0aWMgdm9pZCBtYXRyaXhfZGVhbGxvY2F0ZShNYXRyaXggKm0pCnsKCWZvciAoIHNpemVfdCB5ID0gMDsgeSA8IG0tPnlfZGltOyArK3kgKQoJCWZyZWUobS0+Y2VsbHNbeV0pOwoJZnJlZShtLT5jZWxscyk7CgoJbS0+eF9kaW0gPSAwOwoJbS0+eV9kaW0gPSAwOwoJbS0+Y2VsbHMgPSAwOwp9Cgp2b2lkIG1hdHJpeF9mcmVlKE1hdHJpeCAqbSkKewoJbWF0cml4X2RlYWxsb2NhdGUobSk7CglmcmVlKG0pOwp9Cgp2b2lkIG1hdHJpeF9hc3NpZ24oTWF0cml4ICpkc3QsIE1hdHJpeCBjb25zdCAqc3JjKQp7CgkvLyBub3RlOiBjb3VsZCBvcHRpbWl6ZSB0aGlzIHRvIG5vdCBhbGxvY2F0ZSBpZiBkaW1zIHRoZSBzYW1lCgltYXRyaXhfZGVhbGxvY2F0ZShkc3QpOwoJTWF0cml4ICpuZXcgPSBtYXRyaXhfY29uc3RydWN0KHNyYy0+eV9kaW0sIHNyYy0+eF9kaW0pOwoJZm9yIChzaXplX3QgeSA9IDA7IHkgPCBzcmMtPnlfZGltOyArK3kpCgkJbWVtY3B5KG5ldy0+Y2VsbHNbeV0sIHNyYy0+Y2VsbHNbeV0sIHNyYy0+eF9kaW0gKiBzaXplb2YgbmV3LT5jZWxsc1swXVswXSk7CgoJaWYgKCAhbmV3ICkgZXhpdChFWElUX0ZBSUxVUkUpOwoJKmRzdCA9ICpuZXc7CglmcmVlKG5ldyk7Cn0KCk1hdHJpeCAqbWF0cml4X211bHRpcGx5KE1hdHJpeCBjb25zdCAqbTEsIE1hdHJpeCBjb25zdCAqbTIpCnsKCWlmICggbTEtPnhfZGltICE9IG0yLT55X2RpbSApCgkJcmV0dXJuIE5VTEw7CgoJTWF0cml4ICpuZXcgPSBtYXRyaXhfY29uc3RydWN0KG0xLT55X2RpbSwgbTItPnhfZGltKTsKCglmb3IgKHNpemVfdCBjb2wgPSAwOyBjb2wgPCBtMi0+eF9kaW07ICsrY29sKQkvLyBFYWNoIGNvbHVtbiBvZiBtMgoJCWZvciAoc2l6ZV90IHJvdyA9IDA7IHJvdyA8IG0xLT55X2RpbTsgKytyb3cpCS8vIEVhY2ggcm93IG9mIG0xCgkJewoJCS8vIGRvIGRvdC1wcm9kdWN0IG9mIHRob3NlIHR3bwoJCQlmbG9hdCBzdW0gPSAwOwoKCQkJZm9yIChzaXplX3QgaWkgPSAwOyBpaSA8IG0yLT55X2RpbTsgKytpaSkKCQkJCXN1bSArPSBtMi0+Y2VsbHNbaWldW2NvbF0gKiBtMS0+Y2VsbHNbcm93XVtpaV07CgoJCQluZXctPmNlbGxzW3Jvd11bY29sXSA9IHN1bTsKCQl9CgoJcmV0dXJuIG5ldzsKfQoKc3RhdGljIHZvaWQgbWF0cml4X3ByaW50ZihjaGFyIGNvbnN0ICpwcm9tcHQsIE1hdHJpeCBjb25zdCAqbWF0KQp7CglwcmludGYoIiUyMC4yMHM6IiwgcHJvbXB0KTsKCWZvciAoc2l6ZV90IHkgPSAwOyB5IDwgbWF0LT55X2RpbTsgKyt5KQoJewoJCWlmICh5ID4gMCApIHByaW50ZigiJTIwLjIwcyAiLCAiIik7CgkJcHJpbnRmKCIoICIpOwoJCWZvciAoc2l6ZV90IHggPSAwOyB4IDwgbWF0LT54X2RpbTsgKyt4KQoJCQlwcmludGYoIiUzLjNmICIsIG1hdC0+Y2VsbHNbeV1beF0pOwoJCXByaW50ZigiKVxuIik7Cgl9Cgp9CmludCBtYWluKCkKewoKICAgIGZsb2F0IG1hdDFfYXJyYXlbM11bMl0gPSB7ezAsIDF9LCB7MywgNH0sIHs2LCA3fX07CiAgICBmbG9hdCBtYXQyX2FycmF5WzJdWzNdID0ge3s1LCAxLCAyfSwgezMsIDQsIDV9fTsgCgogICAgTWF0cml4ICptYXQxID0gbWF0cml4X2NvbnN0cnVjdCgzLCAyKTsKICAgIE1hdHJpeCAqbWF0MiA9IG1hdHJpeF9jb25zdHJ1Y3QoMiwgMyk7CiAgICBtYXRyaXhfc2V0X3ZhbHVlcyhtYXQxLCAoZmxvYXQgKikmbWF0MV9hcnJheSk7CiAgICBtYXRyaXhfc2V0X3ZhbHVlcyhtYXQyLCAoZmxvYXQgKikmbWF0Ml9hcnJheSk7CgogICAgTWF0cml4ICptYXQzID0gbWF0cml4X211bHRpcGx5KG1hdDEsIG1hdDIpOwoKCW1hdHJpeF9wcmludGYoIm1hdDEiLCBtYXQxKTsKCW1hdHJpeF9wcmludGYoIm1hdDIiLCBtYXQyKTsKCW1hdHJpeF9wcmludGYoIm1hdDMiLCBtYXQzKTsKCiAgICBtYXRyaXhfb3BlcmF0ZV9zZWxmKG1hdHJpeF9tdWx0aXBseSwgbWF0MSwgbWF0Mik7Cglhc3NlcnQoIG1hdHJpeF9pc19lcXVhbChtYXQxLCBtYXQzKSApOwoKICAgIG1hdHJpeF9mcmVlKG1hdDEpOyAgICAKICAgIG1hdHJpeF9mcmVlKG1hdDIpOyAgICAKICAgIG1hdHJpeF9mcmVlKG1hdDMpOyAgICAKfQo=