#ifndef MATRIX_H
#define MATRIX_H

#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>

#include "common.h"
#include "curves.h"

typedef struct matrix matrix;
struct matrix {
	size_t nrows, ncolumns;
	double * data;
};

static inline double *
rm_pos(const struct matrix * m, size_t r, size_t c)
{
	DEBUG_ASSERT(r < m->nrows);
	DEBUG_ASSERT(c < m->ncolumns);

	return &m->data[r*m->ncolumns+c];
}

static inline double
rm_elem(const struct matrix * m, size_t r, size_t c)
{
	return *rm_pos(m, r, c);
}

static inline double *
cm_pos(const struct matrix * m, size_t r, size_t c)
{
	DEBUG_ASSERT(r < m->nrows);
	DEBUG_ASSERT(c < m->ncolumns);

	return &m->data[c*m->nrows+r];
}

static inline double
cm_elem(const struct matrix * m, size_t r, size_t c)
{
	return *cm_pos(m, r, c);
}

static inline double *
hm_pos(const struct matrix * m, size_t r, size_t c)
{
	DEBUG_ASSERT(r < m->nrows);
	DEBUG_ASSERT(c < m->ncolumns);

	return &m->data[hil_d_from_xy(r,c)];
}

static inline double
hm_elem(const struct matrix * m, size_t r, size_t c)
{
	return *hm_pos(m, r, c);
}

static inline double *
zm_pos(const struct matrix * m, size_t r, size_t c)
{
	DEBUG_ASSERT(r < m->nrows);
	DEBUG_ASSERT(c < m->ncolumns);

	return &m->data[mor_z_from_xy(r,c)];
}

static inline double
zm_elem(const struct matrix * m, size_t r, size_t c)
{
	return *zm_pos(m, r, c);
}

void
matrix_init(struct matrix * m, size_t nrows, size_t ncolumns);

void
matrix_fini(struct matrix * m);

struct matrix *
matrix_create(size_t nrows, size_t ncolumns);

void
matrix_destroy(struct matrix * m);

bool
matrix_equals(const struct matrix * a, const struct matrix * b);

void
matrix_populate(struct matrix * m);

struct matrix *
matrix_convert(const struct matrix * m,
	double (*elem)(const struct matrix * m, size_t r, size_t c),
	double * (*pos)(const struct matrix * m, size_t r, size_t c));

static inline struct matrix *
matrix_rm_to_cm(const struct matrix * m)
{
	return matrix_convert(m, rm_elem, cm_pos);
}

static inline struct matrix *
matrix_rm_to_hm(const struct matrix * m)
{
	return matrix_convert(m, rm_elem, hm_pos);
}

static inline struct matrix *
matrix_rm_to_zm(const struct matrix * m)
{
	return matrix_convert(m, rm_elem, zm_pos);
}

static inline struct matrix *
matrix_cm_to_rm(const struct matrix * m)
{
	return matrix_convert(m, cm_elem, rm_pos);
}

static inline struct matrix *
matrix_cm_to_hm(const struct matrix * m)
{
	return matrix_convert(m, cm_elem, hm_pos);
}

static inline struct matrix *
matrix_cm_to_zm(const struct matrix * m)
{
	return matrix_convert(m, cm_elem, zm_pos);
}

static inline struct matrix *
matrix_hm_to_rm(const struct matrix * m)
{
	return matrix_convert(m, hm_elem, rm_pos);
}

static inline struct matrix *
matrix_hm_to_cm(const struct matrix * m)
{
	return matrix_convert(m, hm_elem, cm_pos);
}

static inline struct matrix *
matrix_hm_to_zm(const struct matrix * m)
{
	return matrix_convert(m, hm_elem, zm_pos);
}

static inline struct matrix *
matrix_zm_to_rm(const struct matrix * m)
{
	return matrix_convert(m, zm_elem, rm_pos);
}

static inline struct matrix *
matrix_zm_to_cm(const struct matrix * m)
{
	return matrix_convert(m, zm_elem, cm_pos);
}

static inline struct matrix *
matrix_zm_to_hm(const struct matrix * m)
{
	return matrix_convert(m, zm_elem, hm_pos);
}

void
matrix_show(const struct matrix * m, FILE * fp);

struct matrix *
matrix_multiply_rm(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_trm(const struct matrix * a, const struct matrix * b, size_t bsize);

struct matrix *
matrix_multiply_cm(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_hm(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_hm2(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_hm3(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_zm(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_atlas(const struct matrix * a, const struct matrix * b);

struct matrix *
matrix_multiply_verify(const struct matrix * a, const struct matrix * b);

#endif