codesamples/C/matrix.c

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include "FEVtypes.h"
#include "FEV.h"


/***********simplified versions with ranges [0,N-1] **************/
/*Vectors are columns*/

dvec dvector(int n){
	return  (dvec) malloc(n*sizeof(double));
}

void dvec_free(dvec v){
	free((char*) v);
}

/* dvec is double-precision vector data, used below for rows. dmat is
double-precision matrix data, the third field of a matrix struct.*/


/* create matrix by allocating r  c-long vectors for the rows,
   a vector of pointers to rows, and a pointer to the latter vector. */

matrix_t mat_new(int r,int c)  /*r rows, c cols*/{
	int i;
	dmat m;
	matrix_t a;
	m = (dmat) malloc( r*sizeof(double *));
	if (m == NULL) return NULL;
	for (i= 0; i<r; i++)
	{
		m[i] = (dvec) malloc( c*sizeof(double));
		if (m[i] == NULL) return NULL;
	}
	a =(matrix_t) malloc ( sizeof(struct matrix));
	if ( a == NULL) return NULL;
	a->r = r;
	a->c = c;
	a->mat = m;
	a-> fev_type = MATRIX;
	return a;
}

/*  can't just free the pointer to the matrix, there's no garbage collection */


void mat_free(matrix_t m){
	int i,r;
	r =m->r;
	for (i = 0; i<r; i++) free ((char *) m->mat[i]);
	free((char *) m->mat);
	free((char *) m);
}

/* not much of an error function...*/
void cberror (char *msg)
{
        fprintf(stderr, "\n %s \n", msg);
}


/* print out a matrix in a human-readable form */
int mat_pr(char *msg, matrix_t m, FILE *fout) {
	int i,j;
	int r,c;
	if (msg == NULL) {
		cberror("null message arg in mat_pr");
		return NILARG1;
	}
	if (m == NULL) {
		cberror("null matrix arg in mat_pr");
		return NILARG2;
	}
	r = m->r;
	c = m->c;
	if(fout == NULL){
		printf("MATRIX: %s Type: %d , Dims %d, %d\n", msg,m->fev_type,r,c);
		for (i=0; i<r; i++){
			for (j=0; j<c; j++)
				printf("%.2f ",m->mat[i][j]);
			printf("\n");
		}
		printf("~\n");
		return OK;
	} else {
		fprintf(fout, "MATRIX: %s Type: %d , Dims %d, %d\n", msg,m->fev_type,r,c);
		for (i=0; i<r; i++){
			for (j=0; j<c; j++)
				fprintf(fout, "%.2f ",m->mat[i][j]);
			fprintf(fout, "\n");
		}
		fprintf(fout, "~\n");
		return OK;
	}
}

/* create a matrix from the output of mat_pr, stored in a file (fname) */
matrix_t mat_frd(char* fname){
	int i, j;
	float n;
	int r, c, assigned;
	matrix_t m;
	FILE *fmat;

	fmat = fopen(fname, "r");
	assigned = fscanf(fmat, "MATRIX: %*s Type: %*d , Dims %d, %d\n",
					 &r, &c);
	if (assigned != 2) {
		cberror("incorrectly formatted matrix file");
		return (matrix_t) NULL;
	}
	m = mat_new(r, c);
	for (i = 0; i < r; i++) {
		for (j=0; j<c; j++) {
			/* get each number */
			assigned = fscanf(fmat, "%f ", &n);
			if(assigned != 1) {
				cberror("Formatting error");
				return (matrix_t) NULL;
			}
			e(m, i, j) = n;
		}
		fscanf(fmat, "\n"); //end of each row
	}
	fclose(fmat);
	return m;

}


/* multiply two matrices a and b; return the result as a new matrix */
matrix_t mat_mpy(matrix_t a, matrix_t b){
	int i, j, k;
	int r, c, n;
	matrix_t m; //output matrix

	/* dims of output matrix*/
	r = a->r;
	c = b->c;

	if (a->c != b->r) {
		cberror("incompatible matrices; could not multiply");
		return (matrix_t) NULL;
	}

	n = a->c;
	m = mat_new(r, c);
	for (i=0; i<r; i++) {
		for (j=0; j<c; j++) {
			for (k=0; k<n; k++) {
				e(m, i, j) += e(a, i, k)*e(b, k, j);
			}
		}
	}
	return m;
}


/* add two matrices a and b: multiply b by scalar sign,
  add to a, and return result */
matrix_t mat_add(matrix_t a, matrix_t b, int sign){
	int i, j;
	int r, c;
	matrix_t m;

	/* dims of output matrix - same as a and b */
	r = a->r;
	c = b->c;
	if ((r != b->r) || (c != a->c)) {
		cberror("matrices of different size; could not add");
		return (matrix_t) NULL;
	}

	m = mat_new(r, c);
	for (i=0; i<r; i++) {
		for (j=0; j<c; j++) {
			m->mat[i][j] = (e(b, i, j)*sign) + e(a, i, j);
		}
	}

	return m;
}