/**********************************************************************
 *                                                                    *
 * This C language library provides some simple matrix operations and *
 * the functions MLOAD and MSAVE enables the use of MAT files         *
 * generated by MATLAB.                                               *
 *                                                                    *
 * Coded by     : Steffen Torp & Magnus Norgaard                      *
 *                                                                    *
 * References   : Brugervejledning til C/Matlab matrix bibliotek,     *
 *                Appendix F i Noergaard og Torp (1993).              *
 *                                                                    *
 * Dependencies : matrix.h                                            *
 *                                                                    *
 * Comments     : All matrices are represented by a structure with    *
 *                the dimensions of the matrix and a pointer to       *
 *                an array of pointers to each row in the matrix.     *
 *                The functions using stdio can be excluded by        *
 *                defining the symbol NO_IO to the compiler.          *
 *                If the library is compiled on the OS/9 system       *
 *                the symbol OS9 must be defined.                     *
 *                                                                    *
 * Last edit    : sept. 22, 1994 Magnus Norgaard                      *
 *                                                                    *
 **********************************************************************/

#include <stdlib.h>
#include <string.h>

#include "matrix2.h"     

#define TRUE 1
#define FALSE 0

#ifndef NO_IO  /* I/O functions are not available in the DSP C compiler */
#include <stdio.h>
int  loadmat(FILE*, int*, char*, int*, int*, int*, double**, double**);
void savemat(FILE*,  int, char*,  int,  int,  int,  double*,  double*);
long longconv(long);           /* Internal functions, only used by  */
void matconv(int, double*);    /* The functions in this library.    */
typedef struct {               /* Matrix structure for MATLAB files */
	long type;
	long mrows;
	long ncols;
	long imagf;
	long namlen;
} Fmatrix;
#endif


#define READ "rb"
#define APND "ab"
#define SYSTEM 0               /* 0 = PC, 1000 = Apollo and OS/9 */
#ifdef OS9
  #undef READ
  #undef APND
  #undef SYSTEM
  #define SYSTEM 1000
  #define READ "r"
  #define APND "a"
#endif


/*
 * MMAKE :
 * -------
 * This function allocates memory for a matrix of the specified size
 * and assigns the specified dimensions to the allocated matrix.
 * Inputs  : rows - number of rows.
 *           cols - number of columns.
 * Outputs : *ptm - pointer to the new matrix structure.
 */
matrix *mmake( int rows, int cols )
{
	matrix *ptm;
	double **row_pointers;
	register int i;

#if RUNCHK

	/* Check that no dimension is zero. */
	if ((rows==0)||(cols==0)) merror("Invalid dimension error in mmake!");
#endif

	/* Memory allocation for matrix structure. */
	ptm = (matrix*)malloc(sizeof(matrix));

	/* Memory for array of pointers to rows. */;
	row_pointers = (double**)malloc(rows*sizeof(double*));

	/* Memory for all rows, initialize row pointers. */
	row_pointers[0] = (double*)malloc(rows*cols*sizeof(double));
	for (i=1;i<rows;i++) {
		row_pointers[i] = row_pointers[i-1] + cols;
	}

#if RUNCHK

	/* Check if last allocation was ok ! */
	if (!row_pointers[0]) {
		merror("Memory allocation error in mmmake!");
	}
#endif

	ptm->row = rows;             /* Initialize matrix structure */
	ptm->col = cols;
	ptm->mat = row_pointers;     /* Pointer to row pointers     */

	return ptm;           /* Return pointer to matrix structure */
}


/*
 * MFREE :
 * -------
 * This function deallocates the memory that was used for a matrix.
 * Input : Pointer to the matrix to be deallocated.
 */
void mfree( matrix *ptm )
{
	/* Deallocate rows */
	free(ptm->mat[0]);

	/* Deallocate row pointer array */
	free(ptm->mat);

	/* Deallocate matrix structure */
	free(ptm);
}



/*
 * MSAVE :
 * -------
 * This function is used to save a matrix in a MAT-file, that can be read
 * by MATLAB. The function uses the SAVEMAT.C function provided by MATLAB.
 * The function only handles numrical matrices.
 * Inputs : ptm       - pointer to the matrix that will be saved.
 *          file_name - string with name of file to save in.
 *          mat_name  - string with name of matrix in MATLAB.
 */
void msave( matrix *ptm, char file_name[], char mat_name[] )
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	FILE *fp;
	double *preal,*pimag;
	double *data;
	int type,imagf,mn,i,j,k;

	mn = ptm->row*ptm->col;                       /* Number of data.  */
	data = (double*)malloc(mn*sizeof(double));    /* Allocate buffer. */
	k = 0;

	/* Prepare for SAVEMAT call using MATLAB format. */

	for ( j=0; j < ptm->col; j++ )
	{
		for ( i=0; i < ptm->row; i++ )
		{
			*(data+k) = get_val(ptm,i,j); /* Write to buffer. */
			k++;
		}
	}
	imagf = 0;                  /* No complex data.                   */
	type  = SYSTEM;             /* Only numeric data.                 */
	preal = data;               /* Real data is in buffer "data".     */
	pimag = (double *)0;        /* No complex data.                   */
	fp = fopen(file_name,APND); /* Append to existing file or create. */
				    /* new one, using binary file-format. */ 
	if (!(fp)) merror("Can't open file error in msave!");
	savemat(fp,type,mat_name,ptm->row,ptm->col,imagf,preal,pimag);
				    /* Save the matrix in MAT-file.       */
	fclose(fp);                 /* Close MAT-file.                    */
	free(data);                 /* Deallocate buffer memory.          */
#else
	ptm=ptm;
	file_name=file_name;
	mat_name=mat_name;
	exit(0); /* printf is not available in DSP compiler. */
#endif
}


/*
 * MLOAD :
 * -------
 * This function is used to load a matrix from a MAT-file, that was written
 * by MATLAB. The function uses the LOADMAT.C function provided by MATLAB,
 * but it only handles numerical and real matrices.
 * Inputs : file_name - string with the name of the MAT-file to read.
 *          mat_name  - string with the name of the matrix to load.
 * Output : ptm       - Pointer to the loaded matrix.
 */
matrix *mload( char file_name[], char mat_name[] )
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	FILE *fp;
	matrix *ptm;
	char name[20];
	int found,i,j,k;
	int type,mrows,ncols,imagf;
	double *xr,*xi;

	fp = fopen(file_name,READ);        /* Open MAT-file for reading   */
					   /* from the beginning of file. */
	if (!(fp)) merror("File not found error in mload!");

	do  /* Search for the matrix with name in mat_name */
	{
		/* Read the file until matrix is found or EOF is reached. */
		if (loadmat(fp,&type,name,&mrows,&ncols,&imagf,&xr,&xi)) {
			printf("The searched matrix is : %s \n",mat_name);
			merror("matrix not in file error in mload!");
		}

		found = strncmp(mat_name,name,20);
	}
	while ( found!=0 );      /* Keep searching until matrix is found. */

	if (imagf) {
		printf("%s\n","Warning from mload :");
		printf("%s\n","The matrix was complex, but the");
		printf("%s\n","imaginary part has been ignored !");
	}

	if ((mrows==0)||(ncols==0)) merror("Matrix was empty error in mload!");

	ptm = mmake( mrows, ncols ); /* Allocate memory for loaded matrix */

	/* Copy data from buffer to matrix using C matrix structure. */
	k = 0;
	for ( j=0; j<ncols; j++ )
	{
		for ( i=0; i<mrows; i++ )
		{                               
			put_val(ptm,i,j,*(xr+k));
			k++;
		}
	}

	free(xr);   /* Deallocate data buffers and pointers.*/
	if (imagf) free(xi);
	fclose(fp); /* Close MAT-file.                      */
	return ptm; /* Return pointer to the loaded matrix. */
#else
	file_name=file_name;
	mat_name=mat_name;
	exit(0); /* printf is not available in DSP compiler. */
#endif
}


/*
 * SLOAD :
 * -------
 * This function loads a scalar (one by one matrix) from a MAT-file.
 * Inputs : file_name - Name of MAT-file to search for scalar.
 *          scal_name - Name of Matlab variable in MAT-file.
 * Output : scalar    - The scalar, which was loaded.
 */
double sload( char file_name[], char scal_name[] )
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	matrix *ptm;
	double scalar;

	ptm    = mload( file_name, scal_name );
	scalar = ptm->mat[0][0];
	mfree(ptm);

	return scalar;
#else
	file_name=file_name;
	scal_name=scal_name;
	exit(0); /* printf is not available in DSP compiler. */
#endif
}


/*
 * MERROR :
 * --------
 * This function displays the specified error message and exits to server
 * operating system.
 * Input : Text string with the error message.
 */
void merror( char error_text[] )
{

#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	printf("--- Runtime error in matrix library:\n");
	printf("--- ");
	printf("%s\n\n",error_text);
#else
	error_text=error_text; /* Dummy command to avoid warnings on DSP */ 
#endif
	/*exit(1);*/
	mexErrMsgTxt("AN ERROR OCCURED IN A CMEX-PROGRAM");
}


/*
 * MPRINT :
 * --------
 * This is a simple function to print a matrix on the screen.
 * Input : ptm - Pointer to the matrix that will be printed.
 */
void mprint( matrix *ptm )
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	int i,j;
	for ( i=0; i < ptm->row; i++ )
	{
		for( j=0; j < ptm->col; j++ )
		{
			printf("%4f %s",ptm->mat[i][j]," ");
		}
		printf("\n");
	}
#else
	ptm=ptm;
	exit(0); /* printf is not available in DSP compiler. */
#endif
}


/*
 * MINPUT1 :
 * ---------
 * This function reads the dimensions and content of a matrix from 
 * the keyboard.
 * Output : Pointer to new matrix.
 */
matrix *minput1()
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	matrix *ptm;
	int i,j,mrows,ncols;

	printf("%s","Number of rows    :");
	scanf("%d",&mrows);
	printf("%s","Number of columns :");
	scanf("%d",&ncols);

	ptm = mmake(mrows,ncols);

	for ( i=0; i<mrows; i++ )
	{
		for ( j=0; j<ncols; j++ )
		{
			printf("%s %d %s %d %s","Enter element (",i,",",j,")  :");
			scanf("%lf",&(ptm->mat[i][j]));
			/* PS. scanf does not work with floating formats in Turbo C. */
		}
	}
	return ptm;
#else
	exit(0); /* scanf is not available in DSP compiler. */
#endif
}


/*
 * MINPUT2 :
 * ---------
 * This function prompts for the content of a matrix of known size.
 * Input : Pointer to matrix to be read from keyboard.
 */
void minput2( matrix *ptm )
{
#ifndef NO_IO     /* For compilation by the DSP ANSI-C compiler. */

	int i,j;

	printf("%s %d %s %d %s\n","This matrix has",ptm->row,"rows and",
			ptm->col,"columns.");
	printf("%s %d %s\n","Please enter the following",
			ptm->row*ptm->col,"elements.");

	for ( i=0; i < ptm->row; i++ )
	{
		for ( j=0; j < ptm->col; j++ )
		{
			printf("%s %d %s %d %s","Enter element (",i,",",j,")  :");
			scanf("%lf",&(*ptm).mat[i][j]);
			/* PS. scanf does not work with floating formats in Turbo C. */
		}
	}
#else
	ptm=ptm;
	exit(0); /* scanf is not available in DSP compiler. */
#endif
}


/*
 * MTRANS :
 * --------
 * This function computes the transposed of a matrix.
 * Inputs : ptm2 - Pointer to argument matrix.
 *          ptm1 - Pointer to result matrix.
 * Notice : The two matrices may NOT be the same.
 */
void mtrans( matrix *ptm1, matrix *ptm2 )
{
	register int i,j;

#if RUNCHK

	if (ptm1==ptm2) merror("Memory conflict error in mtrans!");
#endif

	for ( i=0; i < ptm2->col; i++ )          /* Transpose matrix */
	{
		for ( j=0; j < ptm2->row; j++ )
		{
			ptm1->mat[i][j] = ptm2->mat[j][i];
		}
	}
}


/*
 * MINIT :
 * -------
 * This function initializes a matrix to all zeros.
 * Inputs : ptm - Pointer to the matrix that will be initialized.
 */
void minit( matrix *ptm )
{
	register unsigned char zero_byte = 0;
	register int bytes;
	
	/* Number of bytes in matrix data area. */
	bytes = 8*(ptm->row)*(ptm->col);

	/* Set data area to all zeroes. */
	memset( ptm->mat[0], zero_byte, bytes );
}



/*
 * MINITX :
 * --------
 * This function initializes a matrix to all x's. x being a scalar.
 * Inputs : ptm - Pointer to the matrix that will be initialized.
 *          x   - Element to be put into the matrix  
 */
void minitx( matrix *ptm, double x )
{
	register int i,j;
        for(i=0; i < ptm->row; i++)
	{
	  for(j=0; j < ptm->col; j++)
          {	
	    ptm->mat[i][j]=x;
	  }
	}
}


/*
 * MRAND :
 * -------
 * This function initializes a matrix with random numbers between -1 and 1.
 * Inputs : ptm - Pointer to the matrix that will be initialized.
 */
void mrand( matrix *ptm )
{
	register int i,j;
        for(i=0; i < ptm->row; i++){
	  for(j=0; j < ptm->col; j++)
	              ptm->mat[i][j]=rand()*2.0/RAND_MAX - 1.0;
	}
}



/*
 * MFIND :
 * -------
 * This function finds a given element in a matrix and returns
 * the row and coloumn coordinates to the locations. If no
 * instances of the element is found, a zero by zero matrix is
 * returned.
 * 
 * Inputs : ptm    - Pointer to matrix.
 *          elem   - Element to be searched for (double)
 * Output : retmat - 2 coloumn matrix. 1st coloumn contains row indices,
 *                   2nd coloumn contains the colomn indices.
 */
matrix *mfind( matrix *ptm, double elem )
{
        register int i,j,k;
	matrix *tmp, *retmat;
	k=0;
	tmp = mmake((ptm->row)*(ptm->col),2);

	/* Seach entire matrix */
	for(i=0; i < ptm->row; i++)
	  {
	    for(j=0; j < ptm->col; j++)
	      {
		if(ptm->mat[i][j] == elem)
		  {
		    tmp->mat[k][0]=i;
		    tmp->mat[k][1]=j;
                    k++;
		  }
	      }
	  }
	if(k==0)   /* The element is'nt recognized */
	{
	  retmat = mmake(1,1);
	  retmat->row=0;
	  retmat->col=0;
	}          /* Copy coordinates to a [ 2 | # of instances] matrix */
	else
	  {
	  retmat = mmake(k,2);
	  for(i=0; i < k; i++)
	  {
	    retmat->mat[i][0] = tmp->mat[i][0];
	    retmat->mat[i][1] = tmp->mat[i][1];
	  }
	}
	mfree(tmp);
	return retmat;
}



/*
 * MNOFIND :
 * ---------
 * This function finds all entries in a matrix not equal to a given element
 * and returns the row and coloumn coordinates to the locations. If all entries
 * equals the element, a zero by zero matrix is returned.
 * 
 * Inputs : ptm    - Pointer to matrix.
 *          elem   - Element to be searched for (double)
 * Output : retmat - 2 coloumn matrix. 1st coloumn contains row indices,
 *                   2nd coloumn contains the colomn indices.
 */
matrix *mnofind( matrix *ptm, double elem )
{
        register int i,j,k;
	matrix *tmp, *retmat;
	k=0;
	tmp = mmake((ptm->row)*(ptm->col),2);

	/* Seach entire matrix */
	for(i=0; i < ptm->row; i++)
	  {
	    for(j=0; j < ptm->col; j++)