| 1 |
673 |
markom |
/***********************************************************************
|
| 2 |
|
|
|
| 3 |
|
|
matrix.h - matrix operations
|
| 4 |
|
|
|
| 5 |
|
|
Copyright (C) 1991 Dean Rubine
|
| 6 |
|
|
|
| 7 |
|
|
This program is free software; you can redistribute it and/or modify
|
| 8 |
|
|
it under the terms of the GNU General Public License. See ../COPYING for
|
| 9 |
|
|
the full agreement.
|
| 10 |
|
|
|
| 11 |
|
|
**********************************************************************/
|
| 12 |
|
|
/*
|
| 13 |
|
|
|
| 14 |
|
|
Simple matrix operations
|
| 15 |
|
|
Why I am writing this stuff over is beyond me
|
| 16 |
|
|
|
| 17 |
|
|
*/
|
| 18 |
|
|
|
| 19 |
|
|
/*
|
| 20 |
|
|
|
| 21 |
|
|
This package provides the Matrix and Vector data types
|
| 22 |
|
|
|
| 23 |
|
|
The difference between this matrix package and others is that:
|
| 24 |
|
|
Vectors may be accessed as 1d arrays
|
| 25 |
|
|
Matrices may still be accessed like two dimensional arrays
|
| 26 |
|
|
This is accomplished by putting a structure containing the bounds
|
| 27 |
|
|
of the matrix before the pointer to the (array of) doubles (in the
|
| 28 |
|
|
case of a Vector) or before the pointer to an (array of) pointers
|
| 29 |
|
|
to doubles (in the case of a Matrix).
|
| 30 |
|
|
|
| 31 |
|
|
|
| 32 |
|
|
Vectors and matrices are collectively called "arrays" herein.
|
| 33 |
|
|
*/
|
| 34 |
|
|
|
| 35 |
|
|
#define HEADER(a) ( ((struct array_header *) a) - 1 )
|
| 36 |
|
|
|
| 37 |
|
|
#define NDIMS(a) (int)(HEADER(a)->ndims)
|
| 38 |
|
|
#define NROWS(a) (int)(HEADER(a)->nrows)
|
| 39 |
|
|
#define NCOLS(a) (int)(HEADER(a)->ncols)
|
| 40 |
|
|
#define ISVECTOR(a) (NDIMS(a) == 1)
|
| 41 |
|
|
#define ISMATRIX(a) (NDIMS(a) == 2)
|
| 42 |
|
|
|
| 43 |
|
|
/* Note: this structure is prepended at the beginning of a Vector, and causes
|
| 44 |
|
|
the Vector data type to be 32-byte aligned, but not 64-byte aligned.
|
| 45 |
|
|
If this were a problem, filler could be filler[5] (or more) instead.
|
| 46 |
|
|
--Sharon Perl, 12/17/98. */
|
| 47 |
|
|
|
| 48 |
|
|
struct array_header {
|
| 49 |
|
|
unsigned char ndims; /* 1 = vector, 2 = matrix */
|
| 50 |
|
|
unsigned char nrows;
|
| 51 |
|
|
unsigned char ncols;
|
| 52 |
|
|
unsigned char filler;
|
| 53 |
|
|
};
|
| 54 |
|
|
|
| 55 |
|
|
typedef double **Matrix;
|
| 56 |
|
|
typedef double *Vector;
|
| 57 |
|
|
|
| 58 |
|
|
Vector NewVector(); /* int r; (number of rows) */
|
| 59 |
|
|
Matrix NewMatrix(); /* int r, c; (number of rows, number of columns) */
|
| 60 |
|
|
void FreeVector(); /* Vector v; */
|
| 61 |
|
|
void FreeMatrix(); /* Matrix m; */
|
| 62 |
|
|
void PrintVector(); /* Vector v; char *fmt; any a1,a2,a3,a4,a5,a6,a7,a8 */
|
| 63 |
|
|
void PrintMatrix(); /* Matrix m; char *fmt; any a1,a2,a3,a4,a5,a6,a7,a8 */
|
| 64 |
|
|
double InnerProduct(); /* Vector v1, v2 */
|
| 65 |
|
|
void MatrixMultiply(); /* Matrix m1, m2, prod; */
|
| 66 |
|
|
void VectorTimesMatrix(); /* Vector v; Matrix m; Vector prod; */
|
| 67 |
|
|
void ScalarTimesVector(); /* double s; Vector v; Vector prod; */
|
| 68 |
|
|
double QuadraticForm(); /* Vector v; Matrix m; (computes v'mv) */
|
| 69 |
|
|
double InvertMatrix(); /* Matrix input_matrix, result_matrix (returns det) */
|
| 70 |
|
|
Vector SliceVector(); /* Vector v; BitVector rowmask */
|
| 71 |
|
|
Matrix SliceMatrix(); /* Matrix m; Bitvector rowmask, colmask; */
|
| 72 |
|
|
Vector VectorCopy(); /* Vector v; */
|
| 73 |
|
|
Matrix MatrixCopy(); /* Matrix m; */
|
| 74 |
|
|
Vector InputVector(); /* FILE *f; */
|
| 75 |
|
|
Matrix InputMatrix(); /* FILE *f; */
|
| 76 |
|
|
|
| 77 |
|
|
double InvertSingularMatrix(); /* Matrix input, result (returns det) */
|
| 78 |
|
|
Matrix DeSliceMatrix(); /* Matrix m, double fill, BitVector rowmask, colmask;
|
| 79 |
|
|
Matrix result */
|
| 80 |
|
|
void OutputVector();
|
| 81 |
|
|
void OutputMatrix();
|
| 82 |
|
|
void ZeroVector();
|
| 83 |
|
|
void ZeroMatrix(); /* Matrix m; */
|
| 84 |
|
|
void FillMatrix(); /* Matrix m; double fill; */
|
