#ifndef _MATRIX_H
#define _MATRIX_H

#include "general.h"

typedef short int mtype;

/* The following structure useful for matrices in which most elements are 0. */
typedef struct MatrElem { 
  mtype value; /* the value of a (nonzero) element */
  int row, column; /* position of element in the matrix */
  struct MatrElem *next[4]; 
  /* these are pointers to the next (nonzero) elements of the matrix in 
     the four directions: down, up, right, left */
} MatrElem;

/* the following values must not be changed unless the functions performing
   the basic operations on MatrElem are also updated */
#define mDOWN 1
#define mUP 3
#define mRIGHT 0
#define mLEFT 2

typedef struct matrix {
  int nc; /* number of columns */
  int nr; /* number of rows */
  unsigned int type:3; /* indicates the type of matrix representation */
  unsigned int trans:1; /* indicates whether the matrix is transposed */
  /* when the matrix is transposed, the element (i,j) of the matrix is 
     interpreted as the element (j,i). */
  union {
    mtype *elem; /* matrix representation for type = 1 */ 
    /* nr*nc elements in this order: row 1, row 2, row 3, ... */
    /* If trans = 0, it is not transposed: elem(i,j) = elem[i*nc + j].
       Otherwise,    it is     transposed: elem(i,j) = elem[j*nc + i]. */

    mtype **ar; /* matrix representation for type = 2 */

    struct { /* matrix representation for type = 3 */
      MatrElem** row;  
      MatrElem** column;
      /* rows and columns organized as chains; for instance, row[i] contains 
         the address of the first element of the chain on row i. */
    };
  };
} matrix;

/* Avoid modifying directly a matrix object.
   Use instead the following functions:
   - AllocateMatrix (for memory allocation and initialization).
   - GetMatrixEl (to obtain the value of the element (i,j) of the matrix).
   - SetMatrixEl (to write a number to the element (i,j) of the matrix).
   - DeallocateMatrix (to deallocate the memory used by the matrix).
*/

void AllocateMatrix(matrix* m, int nr, int nc);
/* This is the first step when defining a matrix. Specify the number of 
   rows (nr) and the number of columns (nc). Then, the function will allocate
   memory to represent the matrix. Moreover, the data will be initialized with
   zero (by means of calloc).

   The function selects automatically the matrix type and then calls 
   AllocateMatrixType */

mtype GetMatrixEl(const matrix* m, int i, int j); 
/* Returns the element (i,j) of the matrix */

void SetMatrixEl(matrix* m, int i, int j, mtype x);
/* Writes x to the element (i,j) of the matrix */

#define NumberOfRows(m) ((m).trans?(m).nc:(m).nr)

#define NumberOfColumns(m) ((m).trans?(m).nr:(m).nc)

#define MatrixLength(m) ((m).nr*(m).nc)

void DeallocateMatrix(matrix* m);

void CopyMatrix(const matrix* src, matrix* dst);
/* The function assumes that dst has been already allocated ! */

void CopyBlock(int nx, int ny, const matrix* src, int xs, int ys,\
               matrix* dst, int xd, int yd);
/* Copies a block of nx by ny elements of src to dst. Specifically, it 
   copies the elements (xs+i, ys+j) of src to (xd+i, yd+j) of dst for all
   i = 0...nx-1 and j = 0...ny-1. */

void Copy2ZeroColumn(const matrix* src, int i, matrix* dst, int k, int offset);
  /* This function is for type 3 matrices. 
     Copies the allocated elements of column i of src to the column k of dst. 
     Namely, it does dst[offset + j][k] = src[j][i] for all j such that 
     src[j][i] is allocated. Now, IF the column k of dst is zero, this means
     that we copy the column i of src to the column k of dst. 

     Example:  When a Petri net object pn1 is copied to a new Petri net object
     pn, we could copy the column i of the input matrix of pn1 as follows: 
     Copy2ZeroColumn(&(pn1.in),i,&(pn.in),i,0);
  */

void MakeZeroColumn(matrix* dst, int i);
/* Makes zero the column i of dst. */

void MakeZeroRow(matrix* dst, int i);
/* Makes zero the row i of dst. */

void AllocateMatrixType(int type, matrix* m, int nr, int nc);
/* Allocates the memory required for an nr by nc matrix. Does not deallocate
   memory if already allocated. */

MatrElem* NextEl(const matrix* m, int dir, int pos, MatrElem* current);
/* Determines the next allocated element in the direction 'dir' at the row or
   column number 'pos'. The starting point is the element 'current'. */ 
/* If called within a loop, the function will give the right result if 
   the matrix chains are not modified between calls. */
/* 'dir' is one of mRIGHT, mLEFT, mUP, mDOWN; pos is the row or column number*/


matrix array2matrix(int** ar, int nr, int nc);
/* Converts array to matrix */

matrix vector2matrix(int* v, int nr, int nc);
/* Converts vector to matrix */

int** vector2array(int* v, int nr, int nc);
/* Converts vector to array */

int** matrix2array(matrix m); 
/* Converts a matrix object to an array of integers: result[i][j] = elem(i,j)*/

void free2Darray(int** ar, int nr);
/* Use to deallocate an array created with matrix2array or vector2array; nr is the number of rows */

int findIndices(void* vect, int n, int esize, int** answer );
/* Finds the indices of the nonzero elements of 'vect'. Here, 'vect' is a 
   vector of 'n' elements. Each element of 'vect' has the size 'esize'. The
   return value is the number of nonzero indices. Further, the indices are
   stored in 'answer'. Note that *answer is a pointer to memory allocated 
   in this function. Memory is always allocated, even if there are no indices. 

   This function is similar to the find function of Matlab. */


// Changes in the following can affect TestCondition and TestExpression !
#define GREATER 2
#define GREATER_OR_EQUAL 1
#define EQUAL 0
#define LESS_OR_EQUAL -1
#define LESS -2

inline int TestCondition(int v, int type, int w);
/* Checks whether 'v' and 'w' (in this order) are in the relationship 
   described by 'type', where 'type' may have one of the values in the 
   range GREATER, ..., LESS. If the relationship is satisfied, the 
   function returns 1. Else, it returns 0. */

int TestExpression(const int* v, int c, const int* w, int len);
/*  The vectors 'v' and 'w' are compared for the first 'len' elements.
    'c' may have one of the values in the range GREATER, ..., LESS above. 
    The return value is as follows:
    * If 'len' <= 0 or v = 0 or w = 0 the function returns a nonzero value;
    * Else if 'c' is none of the above, the function returns zero;
    * Else if the expression is satisfied for each pair of elements v[i] and 
      w[i], the function returns a nonzero value;
    * Else the function returns zero. */

inline void InvertSignOfMatrix(matrix* m); /* m = -m */

inline void InvertSignOfVector(int* v, int len); /* v = -v */


// matrixmath.h and extendedmatrix.h of spnbox have additional matrix functions


#endif
