/***********************************************************************
   (C) Copyright, 1990 by Dean Rubine, Carnegie Mellon University
    Permission to use this code for noncommercial purposes is hereby granted.
    Permission to copy and distribute this code is hereby granted provided
    this copyright notice is retained.  All other rights reserved.
 **********************************************************************/

/*
 Simple matrix operations
*/

/*

This package provides the Matrix and Vector data types

The difference between this matrix package and others is that:
	Vectors may be accessed as 1d arrays
	Matrices may still be accessed like two dimensional arrays
This is accomplished by putting a structure containing the bounds
of the matrix before the pointer to the (array of) FLOATs (in the
case of a Vector) or before the pointer to an (array of) pointers
to FLOATs (in the case of a Matrix).


Vectors and matrices are collectively called "arrays" herein.
*/

#define HEADER(a)	( ((struct array_header *) a) - 1 )

#define	NDIMS(a)	(HEADER(a)->ndims)
#define NROWS(a)	(HEADER(a)->nrows)
#define NCOLS(a)	(HEADER(a)->ncols)
#define	ISVECTOR(a)	(NDIMS(a) == 1)
#define	ISMATRIX(a)	(NDIMS(a) == 2)

struct array_header {
	unsigned char	ndims;	/* 1 = vector, 2 = matrix */
	unsigned char	nrows;
	unsigned char	ncols;
	unsigned char	filler;
};

typedef FLOAT **Matrix;
typedef FLOAT *Vector;

Vector	NewVector (int);      /* (number of rows) */
Matrix	NewMatrix (int, int); /* (number of rows, number of columns) */
void	FreeVector (Vector);
void	FreeMatrix (Matrix);
FLOAT	InnerProduct (Vector, Vector);
void	MatrixMultiply (Matrix, Matrix, Matrix); /* m1 * m2 = prod */
void	VectorTimesMatrix (Vector, Matrix, Vector); /* v*m = prod */
void	ScalarTimesVector (FLOAT, Vector, Vector); /* s * v = prod */
FLOAT	QuadraticForm (Vector, Matrix); /* computes v'mv */
FLOAT	InvertMatrix (Matrix, Matrix); /* m2 = 1/m1 (returns det) */
Vector	SliceVector (Vector, BitVector);/* v, rowmask */
Matrix	SliceMatrix (Matrix, BitVector, BitVector);/* m, rowmask, colmask */
Vector	VectorCopy (Vector);
Matrix	MatrixCopy (Matrix);


void	PrintVector ();
		     /* Vector, char*, int a1,
		     int a2,int a3,int a4,int a5,int a6,int a7,int a8 */
void	PrintMatrix ();
		     /* Matrix, char*, int a1,
		     int a2,int a3,int a4,int a5,int a6,int a7,int a8 */
			 
#ifndef MAC
void    OutputVector(FILE*, Vector);
void    OutputMatrix(FILE*, Matrix);
Vector	InputVector (FILE*);
Matrix	InputMatrix (FILE*);
#endif

FLOAT	InvertSingularMatrix (Matrix, Matrix);/* input, result (returns det) */
Matrix  DeSliceMatrix (Matrix, FLOAT, BitVector, BitVector, Matrix); 
                      /* m, fill, rowmask, colmask, result */
void    ZeroVector (Vector);
void	ZeroMatrix (Matrix);
void	FillMatrix (Matrix, FLOAT);