--- /dev/null
+/* matrix.c -- basic matrix operations
+ *
+ * author: hackbard@hackdaworld.dyndns.org
+ *
+ */
+
+#include "matrix.h"
+
+int matrix_init(t_matrix *matrix,int outfd) {
+
+ dprintf(outfd,"[matrix] initializing matrix api ...\n");
+
+ matrix->outfd=outfd;
+ matrix->type=0;
+ matrix->m=0;
+ matrix->n=0;
+ matrix->z=0;
+ matrix->m1=NULL;
+ matrix->m2=NULL;
+
+ return M_SUCCESS;
+}
+
+int matrix_alloc_data(t_matrix *matrix) {
+
+ /* there is always one matrix */
+ if((matrix->m==0)||(matrix->n==0)) {
+ dprintf(matrix->outfd,"[matrix] wrong dimensions ...\n");
+ return M_DIM_FAILURE;
+ }
+
+ matrix->m1=malloc(matrix->m*matrix->n*sizeof(double));
+ if(matrix->m1==NULL) {
+ dprintf(matrix->outfd,"[matrix] malloc failed\n");
+ return M_ALLOC_FAIL;
+ }
+
+ /* but often there are more ... */
+ if(matrix->type&MULT) {
+ if(matrix->z==0) {
+ dprintf(matrix->outfd,"[matrix] wrong dimensions ...\n");
+ return M_DIM_FAILURE;
+ }
+
+ matrix->m2=malloc(matrix->n*matrix->z*sizeof(double));
+ if(matrix->m2==NULL) {
+ dprintf(matrix->outfd,"[matrix] malloc failed\n");
+ return M_ALLOC_FAIL;
+ }
+ matrix->res=malloc(matrix->m*matrix->z*sizeof(double));
+ if(matrix->res==NULL) {
+ dprintf(matrix->outfd,"[matrix] malloc failed\n");
+ return M_ALLOC_FAIL;
+ }
+ } else {
+ dprintf(matrix->outfd,"[matrix] type of operation not supported by now\n");
+ return M_NOT_SUPPORTED;
+ }
+
+ return M_SUCCESS;
+}
+
+int matrix_shutdown(t_matrix *matrix) {
+
+ dprintf(matrix->outfd,"[matrix] shutdown\n");
+
+ if(matrix->m1!=NULL) free(matrix->m1);
+ if(matrix->m2!=NULL) free(matrix->m2);
+ if(matrix->res!=NULL) free(matrix->res);
+
+ return M_SUCCESS;
+}
+
+/* matrix operations */
+
+int matrix_calc(t_matrix *matrix) {
+
+ unsigned int i,j,k,off,off1,off2;
+
+ if(!(matrix->type&MULT)) {
+ dprintf(matrix->outfd,"[matrix] only MULT supported by now ...\n");
+ return M_NOT_SUPPORTED;
+ }
+
+ /* mult operation, an O(N^3) problem
+ *
+ * maybe speedup by different order of mult?
+ *
+ */
+ for(i=0;i<matrix->z;i++) {
+ off1=0;
+ off=i;
+ for(j=0;j<matrix->n;j++) {
+ off+=matrix->z;
+ off2=i;
+ matrix->res[off]=0;
+ for(k=0;k<matrix->m;k++) {
+ matrix->res[off]+=(matrix->m1[off1]*matrix->m2[off2]);
+ off1+=1;
+ off2+=matrix->z;
+ }
+ }
+ }
+
+ return M_SUCCESS;
+}
--- /dev/null
+/* matrix.h -- matrix headers */
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+/* includes */
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+
+/* defines */
+#define M_SUCCESS 1
+#define M_ERROR -1
+#define M_NOT_SUPPORTED -2
+#define M_DIM_FAILURE -3
+#define M_ALLOC_FAIL -4
+
+typedef struct s_matrix {
+ int outfd;
+ unsigned char type;
+#define MULT (1<<0)
+#define ADD (1<<1)
+#define INV (1<<2)
+#define DIAG (1<<3)
+ int m,n,z;
+ double *m1;
+ double *m2;
+ double *res;
+} t_matrix;
+
+/* function prototypes */
+int matrix_init(t_matrix *matrix,int outfd);
+int matrix_alloc_data(t_matrix *matrix);
+int matrix_calc(t_matrix *matrix);
+int matrix_shutdown(t_matrix *matrix);
+
+#endif