- added backwards calculation of hkl from four circle angles. This

required inclusion of a matrix package.
- modified  counter error handling to send a Stop when the _BAD_BUSY
  error is received.
- added an environment interface to the general controller stuff in choco.*
  Also added setting a parameter directly at the controller object.
- Added a driver for the ETH High Temperature Furnace to be used at
  SANS.
This commit is contained in:
cvs
2000-07-12 12:01:19 +00:00
parent 006f10741c
commit d782d43951
44 changed files with 3199 additions and 25 deletions

20
matrix/DEMO.DAT Normal file
View File

@ -0,0 +1,20 @@
1.1 2 3 4 5 6 7 8 9 9
2 3.1 4 5 6 7 8 9 9 8
3 4 5.1 6 7 8 9 9 8 7
3 4 5 6.1 7 8 9 9 8 7
1 2 3 4 5 7.2 7 8 9 9
1 2 3 4 1 6 7 8.3 9 9
6 2 3 4 5 6 7 8 9.8 9
1 2 3 3 5 6 7 8 9 9.8
1 2 3 4 6 6 7 8 9 9
6 2 3 4 5 4 7 8 9 9
3
1
4
1
2
3
1
1
2
6

BIN
matrix/DEMO.EXE Normal file

Binary file not shown.

62
matrix/MAKEFILE.MSC Normal file
View File

@ -0,0 +1,62 @@
#------------------------------------------------------------------------------
# file: makefile.msc
# desc: makefile for matrix library package under Microsoft C/C++ v7.0
# by: patrick ko
# date: 16 Apr 94
#
# note: a slim matrix library matrix.lib will be generated
#------------------------------------------------------------------------------
OBJS = matcreat.obj matdump.obj materr.obj matadd.obj matsub.obj \
matmul.obj matinv.obj matsolve.obj matdet.obj mattran.obj matdurbn.obj \
mattoepz.obj matsubx.obj
CC = cl -c
C = cl
demo.exe: demo.c matrix.lib
$(C) demo.c matrix.lib
matrix.lib: $(OBJS)
lib matrix +matcreat +matdump +materr +matadd.obj ,,
lib matrix +matsub +matmul +matinv +matsolve +matdet.obj ,,
lib matrix +mattran +matdurbn +mattoepz +matsubx ,,
matcreat.obj: matcreat.c matrix.h
$(CC) matcreat.c
matdump.obj: matdump.c matrix.h
$(CC) matdump.c
materr.obj: materr.c matrix.h
$(CC) materr.c
matadd.obj: matadd.c matrix.h
$(CC) matadd.c
matsub.obj: matsub.c matrix.h
$(CC) matsub.c
matmul.obj: matmul.c matrix.h
$(CC) matmul.c
matinv.obj: matinv.c matrix.h
$(CC) matinv.c
matsolve.obj: matsolve.c matrix.h
$(CC) matsolve.c
mattran.obj: mattran.c matrix.h
$(CC) mattran.c
matdet.obj: matdet.c matrix.h
$(CC) matdet.c
matdurbn.obj: matdurbn.c matrix.h
$(CC) matdurbn.c
mattoepz.obj: mattoepz.c matrix.h
$(CC) mattoepz.c
matsubx.obj: matsubx.c matrix.h
$(CC) matsubx.c

74
matrix/MAKEFILE.TC Normal file
View File

@ -0,0 +1,74 @@
#------------------------------------------------------------------------------
# file: makefile.tc
# desc: makefile for matrix library package under Turbo C v2.0
# by: patrick ko
# date: 16 Apr 94
#
# note: a slim matrix library matrix.lib will be generated
#------------------------------------------------------------------------------
DRIVE = C:
OBJS = matcreat.obj matdump.obj materr.obj matadd.obj matsub.obj \
matmul.obj matinv.obj matsolve.obj matdet.obj mattran.obj matdurbn.obj \
mattoepz.obj matsubx.obj
LIB = -L$(DRIVE)\TC\LIB
INC = -I$(DRIVE)\TC\INCLUDE
CC = tcc -c -mh $(INC)
C = tcc -mh $(LIB) $(INC)
LL = tlib matrix.lib
demo.exe: demo.c $(OBJS)
$(C) demo.c matrix.lib
matcreat.obj: matcreat.c matrix.h
$(CC) matcreat.c
$(LL) -+matcreat.obj
matdump.obj: matdump.c matrix.h
$(CC) matdump.c
$(LL) -+matdump.obj
materr.obj: materr.c matrix.h
$(CC) materr.c
$(LL) -+materr.obj
matadd.obj: matadd.c matrix.h
$(CC) matadd.c
$(LL) -+matadd.obj
matsub.obj: matsub.c matrix.h
$(CC) matsub.c
$(LL) -+matsub.obj
matmul.obj: matmul.c matrix.h
$(CC) matmul.c
$(LL) -+matmul.obj
matinv.obj: matinv.c matrix.h
$(CC) matinv.c
$(LL) -+matinv.obj
matsolve.obj: matsolve.c matrix.h
$(CC) matsolve.c
$(LL) -+matsolve.obj
mattran.obj: mattran.c matrix.h
$(CC) mattran.c
$(LL) -+mattran.obj
matdet.obj: matdet.c matrix.h
$(CC) matdet.c
$(LL) -+matdet.obj
matdurbn.obj: matdurbn.c matrix.h
$(CC) matdurbn.c
$(LL) -+matdurbn.obj
mattoepz.obj: mattoepz.c matrix.h
$(CC) mattoepz.c
$(LL) -+mattoepz.obj
matsubx.obj: matsubx.c matrix.h
$(CC) matsubx.c
$(LL) -+matsubx.obj

54
matrix/MAKEFILE.UX Normal file
View File

@ -0,0 +1,54 @@
#------------------------------------------------------------------------------
# file: makefile.ux
# desc: makefile for matrix library package under Unix
# by: patrick ko
# date: 16 Apr 94
#------------------------------------------------------------------------------
OBJS = matcreat.o matdump.o materr.o matadd.o matsub.o \
matmul.o matinv.o matsolve.o mattran.o matdet.o mattoepz.o matdurbn.o
CC = cc -c
C = cc
CO = -lm -g
demo: demo.c $(OBJS)
$(C) demo.c -o demo $(OBJS) $(CO)
matcreat.o: matcreat.c matrix.h
$(CC) matcreat.c
matdump.o: matdump.c matrix.h
$(CC) matdump.c
materr.o: materr.c matrix.h
$(CC) materr.c
matadd.o: matadd.c matrix.h
$(CC) matadd.c
matsub.o: matsub.c matrix.h
$(CC) matsub.c
matmul.o: matmul.c matrix.h
$(CC) matmul.c
matinv.o: matinv.c matrix.h
$(CC) matinv.c
matsolve.o: matsolve.c matrix.h
$(CC) matsolve.c
mattran.o: mattran.c matrix.h
$(CC) mattran.c
matdet.o: matdet.c matrix.h
$(CC) matdet.c
mattoepz.o: mattoepz.c matrix.h
$(CC) mattoepz.c
matdurbn.o: matdurbn.c matrix.h
$(CC) matdurbn.c
matsubx.o: matsubx.c matrix.h
$(CC) matsubx.c

482
matrix/MATRIX.DOC Normal file
View File

@ -0,0 +1,482 @@
/*
*-----------------------------------------------------------------------------
* file: matrix.doc
* desc: document for matrix toolbox function calls
* by: KO shu pui, patrick
* date: 24 may 92 v0.4
* 23 sep 93 v0.41
* 16 apr 94 v0.42
*-----------------------------------------------------------------------------
*/
===============================================================================
0. INTRODUCTION
===============================================================================
This document only provides you the following information:
0.1 How to create and free a matrix
0.2 Description of each matrix function call
0.3 Some hints to use this toolbox
Remember that this document will NOT describe the data structure and
any technical details of the toolbox - just because this document is
aimed at to be a sort-of "User's Guide" for programmers.
===============================================================================
1. OUR MATRIX OF REFERENCE
===============================================================================
In order to avoid terms confusion, here is our matrix of reference
(matrix A) which is of size m x n.
Column
Row 0 1 2 n-1
0 [ a0,0 a0,1 a0,2 ... a0,n-1 ]
A = 1 [ a1,0 a1,1 a1,2 ... a1,n-1 ]
2 [ a2,0 a2,1 a2,2 ... a2,n-1 ]
[ ... ... ... ... ... ]
m-1 [ am-1,0 am-1,1 am-1,2 ... am-1,n-1 ]
===============================================================================
2. BASIC MATRIX OBJECT OPERATION
===============================================================================
-------------------------------------------------------------------------------
Function : MATRIX mat_creat (int m, int n, int type)
Synopsis : creation of a matrix which can be used by the matrix toolbox
functions; memory is allocated for this object; and some
initialization is performed.
Parameter: m: number of rows
n: number of columns
type: matrix initialization type where
type=
UNDEFINED do not initialize the matrix
ZERO_MATRIX fill zero to all matrix elements
UNIT_MATRIX fill a one to all matrix element ai,j
where i=j
Return Value: the matrix object
Example:
MATRIX A;
/*
* create a 15 x 15 matrix;
* do not initialize the elements
*/
A = mat_creat( 15, 15, UNDEFINED);
/*
* put a 4 in element A(0,14) of matrix A,
* put a 2 in element A(3,5) of matrix A
*/
A[0][14] = 4.0;
A[3][5] = 2.0;
See Also: mat_free(), for Accessing a matrix, see this example.
-------------------------------------------------------------------------------
Function: MATRIX mat_fill ( MATRIX A, int type )
Synopsis: initialize a matrix will a simple type
Parameter: A: the matrix object for which mat_creat() has been called
type: matrix initialization type where
type=
UNDEFINED do not initialize the matrix
ZERO_MATRIX fill zero to all matrix elements
UNIT_MATRIX fill a one to all matrix element ai,j
where i=j
Return Value: MATRIX A
Example:
MATRIX A;
...
/*
* fill the matrix A will zeroes
*/
mat_fill( A, ZERO_MATRIX );
See Also: mat_creat()
-------------------------------------------------------------------------------
Function : int mat_free ( MATRIX A )
Synopsis : free all memory occupied by the matrix object A
Parameter: A: the matrix object for which mat_creat() was called
Return Value: None
Example:
MATRIX A;
A = mat_creat(...);
...
mat_free(A);
Note: since memory may be a very scarce resource in a computer,
as a C programmer you are advised to free up the matrix as
soon as that matrix will not be used any more in future.
See Also: mat_creat()
-------------------------------------------------------------------------------
Function: MatCol ( A )
Synopsis: find out the number of columns of a matrix object A
Parameter: A: the matrix object for which mat_creat() was called
Return Value: number of columns
Example:
int i;
...
i = MatCol(A);
Note: this is a macro
See Also: MatRow()
-------------------------------------------------------------------------------
Function: MatRow ( A )
Synopsis: find out the number of rows of a matrix object A
Parameter: A: the matrix object for which mat_creat() was called
Return Value: number of rows
Example:
int i;
...
i = MatRow(A);
Note: this is a macro
See Also: MatCol()
-------------------------------------------------------------------------------
Function: MATRIX mat_colcopy1 ( MATRIX A, MATRIX B, int j1, int j2 )
Synopsis: copy a column from a matrix A to a column at matrix B
Parameter: A, B: the matrix objects for which mat_creat() was called
column j1 of A is copied to column j2 of B;
Return Value: MATRIX A
Example:
MATRIX A, B;
A = mat_creat( 5, 4, ZERO_MATRIX );
B = mat_creat( 5, 4, UNDEFINED );
/*
* copy column 2 of A to column 0 of B
*/
mat_colcopy1( A, 2, B, 0 );
Note: the sizes of rows of A, B must be the same
See Also: mat_copy()
-------------------------------------------------------------------------------
Function: MATRIX mat_copy ( MATRIX A )
Synopsis: duplicate a matrix
Parameter: A: the matrix object for which mat_creat() was called
Return Value: Another allocated matrix object whose contents are same
as A
Example:
MATRIX A, B;
A = mat_creat( 5, 4, ZERO_MATRIX );
/*
* B is also a 5 x 4 zero matrix now
*/
B = mat_copy( A );
...
mat_free( B );
See Also: mat_colcopy1()
-------------------------------------------------------------------------------
===============================================================================
3. BASIC MATRIX OBJECT INPUT/OUTPUT OPERATION
===============================================================================
-------------------------------------------------------------------------------
Function: int fgetmat (MATRIX A, FILE * fp)
Synopsis: read a matrix from stream
Parameter: A: allocated matrix
fp: a stream pointer obtained from fopen() or predefined
pointer in standard library such as stdin
Return Value: number of matrix elements read
Example:
MATRIX A;
FILE *fp;
A = mat_creat(3, 3, UNDEFINED);
fp = fopen("mymatrix.dat", "r");
fgetmat( A, fp );
See Also: mat_dumpf(), mat_dump(), mat_fdump(), mat_fdumpf()
-------------------------------------------------------------------------------
Function: MATRIX mat_dump (MATRIX A)
MATRIX mat_dumpf (MATRIX A, char * format)
MATRIX mat_fdump (MATRIX A, FILE * fp)
MATRIX mat_fdumpf (MATRIX A, char * format, FILE * fp)
Synopsis: mat_dump: dump a matrix to std output with default format
mat_dumpf: dump a matrix to std output with specified format
mat_fdump: dump a matrix to a file with default format
mat_fdumpf:dump a matrix to a file with specified format
Parameter: A: allocated matrix
format: same as printf() 's format parameter, but can only
be floating point type, such as "%.2f ", etc.
fp: file pointer obtained via fopen()
Return Value: matrix A
Example:
MATRIX A;
A = mat_creat( ... );
...
/*
* dump the matrix to standard output and each element
* is restricted to 2 precision format
*/
mat_dumpf( A, "%.2f ");
See Also: fgetmat()
-------------------------------------------------------------------------------
===============================================================================
4. BASIC MATRIX OBJECT MATHEMATICAL OPERATION
===============================================================================
-------------------------------------------------------------------------------
Function: MATRIX mat_add (MATRIX A, MATRIX B);
MATRIX mat_sub (MATRIX A, MATRIX B);
Synopsis: mat_add: addition of 2 matrices
mat_sub: substraction of 2 matrices
Parameter: A, B: allocated matrices
Return Value: a newly allocated matrix which is the result of the operation
Example:
MATRIX A, B, C;
A = mat_creat( 5, 5, UNIT_MATRIX );
B = mat_creat( 5, 5, UNIT_MATRIX );
C = mat_add( A, B );
mat_dump( C );
mat_free( A );
mat_free( B );
mat_free( C );
Note: A, B must be of the same dimensions
See Also: mat_mul, mat_inv
-------------------------------------------------------------------------------
Function: MATRIX mat_mul (MATRIX A, MATRIX B);
Synopsis: multiplication of 2 matrices
Parameter: A, B: allocated matrix
Return Value: a newly allocated matrix which is the result of the operation
Example:
MATRIX A, B, C;
A = mat_creat( 5, 3, UNIT_MATRIX );
B = mat_creat( 3, 6, UNIT_MATRIX );
/*
* C is now a 5 x 6 matrix
*/
C = mat_add( A, B );
mat_dump( C );
...
Note: the column dimension of A must equal row dimension of B
See Also: mat_add, mat_sub
-------------------------------------------------------------------------------
Function: MATRIX mat_inv (MATRIX A)
Synopsis: find the inverse of a matrix
Parameter: A: allocated square matrix
Return Value: a newly allocated square matrix which is the inverse of A
Example:
MATRIX A, AI;
/*
* A must be a square matrix
*/
A = mat_creat( 7, 7, UNIT_MATRIX );
...
/*
* find the inverse of A
*/
AI = mat_inv( A );
...
See Also: mat_tran, mat_add, mat_sub
-------------------------------------------------------------------------------
Function: MATRIX mat_tran( MATRIX A )
Synopsis: find the transpose of a matrix
Parameter: A: allocated matrix
Return Value: a newly allocated matrix which is the transpose of A
Example:
MATRIX A, T;
A = mat_creat( 3, 5, UNDEFINED );
...
T = mat_tran( A );
...
See Also: mat_inv
-------------------------------------------------------------------------------
===============================================================================
5. DETERMINANT OPERATIONS
===============================================================================
-------------------------------------------------------------------------------
Function: MATRIX mat_submat (MATRIX A, int i, int j)
Synopsis: form a new matrix by deleting a row and a column of A
Parameter: A: allocated matrix
i, j: row and column of A which will not appear in the
resulting matrix
Return Value: a new matrix whose dimensions are 1 less than A
Example:
MATRIX A, B
A = mat_creat( 3, 4, UNDEFINED );
...
B = mat_submat( A, 2, 2 );
/*
* suppose A = [ 1 2 3 4 ]
* [ 3 3 4 5 ]
* [ 6 7 9 9 ]
*
* then B = [ 1 2 4 ]
* [ 3 3 5 ]
*/
Note: Do not be misled -- the contents in A will NOT be changed
See Also: mat_det, mat_cofact, mat_minor
-------------------------------------------------------------------------------
Function: double mat_cofact (MATRIX A, int i, int j)
Synopsis: find out the cofactor of A(i,j)
Parameter: A: allocated square matrix
i,j: row, column position of A for the cofactor
Return Value: the value of cofactor of A(i,j)
Example:
MATRIX A;
A = mat_creat( 5, 5, UNIT_MATRIX );
...
printf( "cofactor of A(0,2) = %f\n", mat_cofact( A, 0, 2 ));
See Also: mat_det, mat_minor
-------------------------------------------------------------------------------
Function: double mat_minor (MATRIX A, int i, int j)
Synopsis: find out the minor of A(i,j)
Parameter: A: allocated square matrix
i,j: row, column position of A for the minor
Return Value: the value of minor of A(i,j)
Example:
MATRIX A;
A = mat_creat( 5, 5, UNIT_MATRIX );
...
printf( "minor of A(0,2) = %f\n", mat_minor( A, 0, 2 ));
See Also: mat_det, mat_cofact
-------------------------------------------------------------------------------
Function: double mat_det (MATRIX A)
Synopsis: find the determinant of a matrix
Parameter: A: allocated square matrix
Return Value: the determinant value of |A|
Example:
MATRIX A;
double det;
A = mat_creat( 5, 5, UNIT_MATRIX );
det = mat_det( A );
See Also: mat_cofact, mat_minor, mat_submat
-------------------------------------------------------------------------------
===============================================================================
6. ADVANCED MATRIX OBJECT MATHEMATICAL OPERATION
===============================================================================
-------------------------------------------------------------------------------
Function: MATRIX mat_lsolve (MATRIX A, MATRIX B)
Synopsis: solve simultaneous linear equation AX = B given A, B
Parameter: A, B: allocated matrix
Return Value: newly allocated matrix X in AX = B
Example:
MATRIX A, B, X;
A = mat_creat( 5, 5, UNDEFINED );
fgetmat( A, stdin );
...
B = mat_creat( 5, 1, UNDEFINED );
fgetmat( B, stdin );
...
X = mat_lsolve( A, B );
mat_dump( X );
Note: number of rows in A and B must be equal
See Also: mat_lsolve_durbin
-------------------------------------------------------------------------------
Function: MATRIX mat_lsolve_durbin (MATRIX A, MATRIX B)
Synopsis: simultaneous linear equations wtih Levinson-Durbin method
Parameter: A, B: allocated matrix where A is an autocorrelated matrix and
B is derived from A
A, B must be the following forms:
| a0 a1 a2 .. an-1 | | x1 | | a1 |
| a1 a0 a1 .. an-2 | | x2 | | a2 |
| a2 a1 a0 .. an-3 | | x3 | = | .. |
| ... | | .. | | .. |
| an-1 an-2 .. .. a0 | | xn | | an |
where A is a symmetric Toeplitz matrix and B
in the above format (related to A)
Return Value: a newly allocated matrix X which is the solution of AX = B
Example: this method only applies to certain A, B only.
e.g.
A = [1 3 9] B = [3]
[3 1 3] [9]
[9 3 1] [12] <- the last one is unrelated to A
See Also: mat_SymToeplz
-------------------------------------------------------------------------------
Function: MATRIX mat_SymToeplz (MATRIX R);
Synopsis: create a symmetric Toeplitz matrix from a
Autocorrelation vector
Parameter: R: allocated matrix of dimension n x 1
Return Value: a newly allocated symmetrical Toeplitz matrix
Example:
e.g.
MATRIX A, R;
R = mat_creat( 3, 1, UNDEFINED );
...
A = mat_SymToeplz( R );
/*
* if
*
* R = [1 3 9] A = [1 3 9]
* [3 1 3]
* [9 3 1]
*
*/
See Also: mat_lsolve_durbin
-------------------------------------------------------------------------------

112
matrix/READ.ME Normal file
View File

@ -0,0 +1,112 @@
-------------------------------------------------------------------------------
Small Matrix Toolbox for C programmers
version 0.42
(Support Unix and DOS)
by Patrick KO Shu-pui
Copyright (c) 1992, 1993, 1994 All Rights Reserved.
-------------------------------------------------------------------------------
ADDRESS TO CONTACT:
fidonet: 6:700/132 BiG Programming Club
[852] 663-0223 19.2 kbps
[852] 663-0236 16.8 kbps
internet: pko@hk.super.net
mailing: Patrick Ko
G.P.O. Box 7468
Hong Kong
-------------------------------------------------------------------------------
MATRX042.ZIP contains
READ .ME - this file
DEMO .C - demo how to use this package
DEMO .DAT - demo data
DEMO .EXE - demo executable for DOS
MAKEFILE.MSC - makefile for Microsoft C/C++ 7.0 on DOS
MAKEFILE.TC - makefile for Turbo C 2.0 on DOS
MAKEFILE.UX - makefile for Unix
MATRIX .DOC - matrix toolbox interface document
MATRIX .H - matrix header file (must include it)
MATADD .C - matrix addition
MATCREAT.C - matrix creation
MATDET .C - find minor, cofactor, determinant
MATDUMP .C - matrix dump
MATERR .C - matrix error handling routine
MATINV .C - matrix inversion
MATMUL .C - matrix multiplication
MATSOLVE.C - linear equations solver
MATSUB .C - matrix substraction
MATSUBX .C - submatrix operation
MATTOEPZ.C - create symmetric Toeplitz matrix
MATDURBN.C - Symmetrix Toeplitz matrix fast solving algorithm
MATTRAN .C - matrix transpose
-------------------------------------------------------------------------------
WHATS NEW in v0.1:
- +, -, *, inverse matrix operations
WHATS NEW in v0.2:
- Linear equation solver
- C-programmer-friendly C sources
WHATS NEW in v0.3:
- better data structure (more Object-Oriented)
- finding minors, cofactors, determinants
- Levinson-Durbin algorithm for symmetric Toeplitz matrix
WHATS NEW in v0.4:
- Revised method for minors, cofactors and determinants
whose time complexity is T(n^3) instead of nearly T(n!).
This is important when you want to find the determinant
of a matrix whose size is over 10 x 10.
- submatrix operator
- matrix formmated dump function
- brief matrix toolbox interface document included
WHATS NEW in v0.41:
- bug fix for unit matrix creation
WHATS NEW in v0.42:
- support Microsoft C/C++ 7.0
HOW TO COMPILE:
1. All Unix environment - make -f makefile.ux
2. DOS (Turbo C v2.0) - make -fmakefile.tc
3. DOS (Microsoft C/C++ v7.0) - nmake -fmakefile.msc
REFERENCES
[1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
John Wiley & Sons, 2nd Ed., 1983. Chap 3.
[2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
John Wiley & Sons, 1978.
[3] Shuzo Saito, Kazuo Nakata, "Fundamentals of Speech Signal
Processing," Academic Press, 1985.
[4] Alfred V.Aho, John E.Hopcroft, Jeffrey D.Ullman, "The Design
and Analysis of Computer Algorithms," 1974.
AUTHOR
All the sources are written by Patrick KO Shu Pui
BiG Programming Club (6:700/132, fidonet)
===============================================================================
AUTHORIZATION NOTICE
This C source package MATRX042.ZIP is FREE for ACADEMIC purpose only.
For COMMERCIAL use, authorization is required from the author.
Please send US$25 for registration.
===============================================================================
DISCLAIMER (I hate this but I have to do that)
You are on your own risk - the author is not responsible for any lost due
to the use of this toolbox.
===============================================================================

92
matrix/demo.c Normal file
View File

@ -0,0 +1,92 @@
/*
*-----------------------------------------------------------------------------
* file: demo.c
* desc: demostrate how to use Patrick's matrix toolbox
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi: 14 may 92 v0.2
* 24 may 92 v0.4
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
* [3] Alfred V.Aho, John E.Hopcroft, Jeffrey D.Ullman, "The Design
* and Analysis of Computer Algorithms," 1974.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include <time.h>
#include "matrix.h"
int main()
{
MATRIX A, B, X, M;
FILE *fp;
double result;
time_t t1, t2;
int tinv, tdet, tmul;
A = mat_creat( 10, 10, UNDEFINED );
B = mat_creat( 10, 1, UNDEFINED );
if ((fp = fopen( "demo.dat", "r" )) == NULL)
{
fprintf( stderr, "file cannot be opened\n" );
exit (0);
}
fgetmat( A, fp );
printf( "|- Matrix A -|\n");
mat_dumpf( A, "%+06.1f " );
t1 = time(&t1);
result = mat_det(A);
t2 = time(&t2);
tdet = t2 - t1;
printf( "\n\nDet(A) = %f\n", result );
printf( "|- Inv A -|\n");
t1 = time(&t1);
X = mat_inv( A );
t2 = time(&t2);
tinv = t2 - t1;
if (X == NULL)
printf( "A is a singular matrix\n" );
else
{
mat_dumpf(X, "%+06.1f ");
printf( "|- A x Inv A -|\n");
t1 = time(&t1);
M = mat_mul( X, A );
t2 = time(&t2);
tmul = t2 - t1;
mat_dumpf( M, "%+06.1f " );
mat_free(M);
mat_free(X);
}
fgetmat( B, fp );
printf( "|- Matrix B -|\n");
mat_dumpf( B, "%+06.1f " );
printf( "|- A x B -|\n");
mat_free(mat_dumpf(mat_mul(A, B), "%+06.1f "));
printf( "time for finding 10 x 10 matrix inversion is less than %d secs\n", tinv );
printf( "time for finding 10 x 10 matrix determinant is less than %d secs\n", tdet );
printf( "time for finding 10 x 10 matrix multiplication is less than %d secs\n", tmul );
mat_free( A );
mat_free( B );
fclose(fp);
}


43
matrix/matadd.c Normal file
View File

@ -0,0 +1,43 @@
/*
*-----------------------------------------------------------------------------
* file: matadd.c
* desc: matrix addition
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_add
* desct: addition of two matrice
* given: A, B = Compatible matrice to be added
* retrn: NULL if malloc() fails
* else allocated matrix of A + B
* comen:
*-----------------------------------------------------------------------------
*/
MATRIX mat_add( A, B )
MATRIX A, B;
{
int i, j;
MATRIX C;
if ((C = mat_creat( MatRow(A), MatCol(A), UNDEFINED )) == NULL)
return (NULL);
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(A); j++)
{
C[i][j] = A[i][j] + B[i][j];
}
return (C);
}

198
matrix/matcreat.c Normal file
View File

@ -0,0 +1,198 @@
/*
*-----------------------------------------------------------------------------
* file: matcreat.c
* desc: matrix mathematics - object creation
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi: 14 may 92 v0.2
* 21 may 92 v0.3
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#ifdef __TURBOC__
#include <alloc.h>
#else
#include <malloc.h>
#endif
#include "matrix.h"
MATRIX _mat_creat( row, col )
int row, col;
{
MATBODY *mat;
int i, j;
if ((mat = (MATBODY *)malloc( sizeof(MATHEAD) + sizeof(double *) * row)) == NULL)
return (mat_error( MAT_MALLOC ));
for (i=0; i<row; i++)
{
if ((*((double **)(&mat->matrix) + i) = (double *)malloc(sizeof(double) * col)) == NULL)
return (mat_error( MAT_MALLOC ));
}
mat->head.row = row;
mat->head.col = col;
return (&(mat->matrix));
}
/*
*-----------------------------------------------------------------------------
* funct: mat_creat
* desct: create a matrix
* given: row, col = dimension, type = which kind of matrix
* retrn: allocated matrix (use mat_free() to free memory)
*-----------------------------------------------------------------------------
*/
MATRIX mat_creat( row, col, type )
int row, col, type;
{
MATRIX A;
if ((A =_mat_creat( row, col )) != NULL)
{
return (mat_fill(A, type));
}
else
return (NULL);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_fill
* desct: form a special matrix
* given: A = matrix, type = which kind of matrix
* retrn: A
*-----------------------------------------------------------------------------
*/
MATRIX mat_fill( A, type )
MATRIX A;
int type;
{
int i, j;
switch (type)
{
case UNDEFINED:
break;
case ZERO_MATRIX:
case UNIT_MATRIX:
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(A); j++)
{
if (type == UNIT_MATRIX)
{
if (i==j)
{
A[i][j] = 1.0;
continue;
}
}
A[i][j] = 0.0;
}
break;
}
return (A);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_free
* desct: free an allocated matrix
* given: A = matrix
* retrn: nothing <actually 0 = NULL A passed, 1 = normal exit>
*-----------------------------------------------------------------------------
*/
int mat_free( A )
MATRIX A;
{
int i;
if (A == NULL)
return (0);
for (i=0; i<MatRow(A); i++)
{
free( A[i] );
}
free( Mathead(A) );
return (1);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_copy
* desct: duplicate a matrix
* given: A = matrice to duplicated
* retrn: C = A
* comen:
*-----------------------------------------------------------------------------
*/
MATRIX mat_copy( A )
MATRIX A;
{
int i, j;
MATRIX C;
if ((C = mat_creat( MatRow(A), MatCol(A), UNDEFINED )) == NULL)
return (NULL);
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(A); j++)
{
C[i][j] = A[i][j];
}
return (C);
}
MATRIX mat_colcopy1( A, B, cola, colb )
MATRIX A, B;
int cola, colb;
{
int i, n;
n = MatRow(A);
for (i=0; i<n; i++)
{
A[i][cola] = B[i][colb];
}
return (A);
}
int fgetmat( A, fp )
MATRIX A;
FILE *fp;
{
int i, j, k=0;
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(A); j++)
{
/*
* to avoid a bug in TC
*/
#ifdef __TURBOC__
{
double temp;
k += fscanf( fp, "%lf", &temp );
A[i][j] = temp;
}
#else
k += fscanf( fp, "%lf", &A[i][j] );
#endif
}
return (k);
}

115
matrix/matdet.c Normal file
View File

@ -0,0 +1,115 @@
/*
*-----------------------------------------------------------------------------
* file: matdet.c
* desc: determinant calculations
* by: ko shu pui, patrick
* date: 21 may 92 v0.3
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
static double signa[2] = {1.0, -1.0};
/*
*-----------------------------------------------------------------------------
* funct: mat_minor
* desct: find minor
* given: A = a square matrix,
* i=row, j=col
* retrn: the minor of Aij
*-----------------------------------------------------------------------------
*/
double mat_minor( A, i, j )
MATRIX A;
int i, j;
{
MATRIX S;
double result;
S = mat_submat(A, i, j);
result = mat_det( S );
mat_free(S);
return (result);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_cofact
* desct: find cofactor
* given: A = a square matrix,
* i=row, j=col
* retrn: the cofactor of Aij
*-----------------------------------------------------------------------------
*/
double mat_cofact( A, i, j )
MATRIX A;
int i, j;
{
double result;
result = signa[(i+j)%2] * A[i][j] * mat_minor(A, i, j);
return (result);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_det
* desct: find determinant
* given: A = matrix
* retrn: the determinant of A
* comen:
*-----------------------------------------------------------------------------
*/
double mat_det( a )
MATRIX a;
{
MATRIX A, P;
int i, j, n;
double result;
n = MatRow(a);
A = mat_copy(a);
P = mat_creat(n, 1, UNDEFINED);
/*
* take a LUP-decomposition
*/
i = mat_lu(A, P);
switch (i)
{
/*
* case for singular matrix
*/
case -1:
result = 0.0;
break;
/*
* normal case: |A| = |L||U||P|
* |L| = 1,
* |U| = multiplication of the diagonal
* |P| = +-1
*/
default:
result = 1.0;
for (j=0; j<MatRow(A); j++)
{
result *= A[(int)P[j][0]][j];
}
result *= signa[i%2];
break;
}
mat_free(A);
mat_free(P);
return (result);
}

78
matrix/matdump.c Normal file
View File

@ -0,0 +1,78 @@
/*
*-----------------------------------------------------------------------------
* file: matdump.c
* desc: matrix mathematics - object dump
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi: 14 may 92 v0.2
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_dump
* desct: dump a matrix
* given: A = matrice to dumped
* retrn: nothing
* comen: matrix a dumped to standard output
*-----------------------------------------------------------------------------
*/
MATRIX mat_dump( A )
MATRIX A;
{
return(mat_fdumpf(A, "%f ", stdout));
}
/*
*-----------------------------------------------------------------------------
* funct: mat_dumpf
* desct: dump a matrix with format string to standard output
* given: A = matrice to dumped
* retrn: nothing
* comen: matrix a dumped to standard output
*-----------------------------------------------------------------------------
*/
MATRIX mat_dumpf( A, s )
MATRIX A;
char *s;
{
return (mat_fdumpf(A, s, stdout));
}
MATRIX mat_fdump( A, fp )
MATRIX A;
FILE *fp;
{
return (mat_fdumpf(A, "%f ", fp));
}
MATRIX mat_fdumpf( A, s, fp )
MATRIX A;
char *s;
FILE *fp;
{
int i, j;
for (i=0; i<MatRow(A); i++)
{
for (j=0; j<MatCol(A); j++)
{
fprintf( fp, s, A[i][j] );
}
fprintf( fp, "\n" );
}
return (A);
}


133
matrix/matdurbn.c Normal file
View File

@ -0,0 +1,133 @@
/*
*-----------------------------------------------------------------------------
* file: matdurbn.c
* desc: Levinson-Durbin algorithm
* by: ko shu pui, patrick
* date:
* revi: 21 may 92 v0.3
* ref:
*
* [1] "Fundementals of Speech Signal Processing," Shuzo Saito,
* Kazuo Nakata, Academic Press, New York, 1985.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_durbin
* desct: Levinson-Durbin algorithm
*
* This function solve the linear eqns Ax = B:
*
* | v0 v1 v2 .. vn-1 | | a1 | | v1 |
* | v1 v0 v1 .. vn-2 | | a2 | | v2 |
* | v2 v1 v0 .. vn-3 | | a3 | = | .. |
* | ... | | .. | | .. |
* | vn-1 vn-2 .. .. v0 | | an | | vn |
*
* where A is a symmetric Toeplitz matrix and B
* in the above format (related to A)
*
* given: R = autocorrelated matrix (v0, v1, ... vn) (dim (n+1) x 1)
* retrn: x (of Ax = B)
*-----------------------------------------------------------------------------
*/
MATRIX mat_durbin( R )
MATRIX R;
{
int i, i1, j, ji, p, n;
MATRIX W, E, K, A, X;
p = MatRow(R) - 1;
W = mat_creat( p+2, 1, UNDEFINED );
E = mat_creat( p+2, 1, UNDEFINED );
K = mat_creat( p+2, 1, UNDEFINED );
A = mat_creat( p+2, p+2, UNDEFINED );
W[0][0] = R[1][0];
E[0][0] = R[0][0];
for (i=1; i<=p; i++)
{
K[i][0] = W[i-1][0] / E[i-1][0];
E[i][0] = E[i-1][0] * (1.0 - K[i][0] * K[i][0]);
A[i][i] = -K[i][0];
i1 = i-1;
if (i1 >= 1)
{
for (j=1; j<=i1; j++)
{
ji = i - j;
A[j][i] = A[j][i1] - K[i][0] * A[ji][i1];
}
}
if (i != p)
{
W[i][0] = R[i+1][0];
for (j=1; j<=i; j++)
W[i][0] += A[j][i] * R[i-j+1][0];
}
}
X = mat_creat( p, 1, UNDEFINED );
for (i=0; i<p; i++)
{
X[i][0] = -A[i+1][p];
}
mat_free( A );
mat_free( W );
mat_free( K );
mat_free( E );
return (X);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_lsolve_durbin
* desct: Solve simultaneous linear eqns using
* Levinson-Durbin algorithm
*
* This function solve the linear eqns Ax = B:
*
* | v0 v1 v2 .. vn-1 | | a1 | | v1 |
* | v1 v0 v1 .. vn-2 | | a2 | | v2 |
* | v2 v1 v0 .. vn-3 | | a3 | = | .. |
* | ... | | .. | | .. |
* | vn-1 vn-2 .. .. v0 | | an | | vn |
*
* domain: where A is a symmetric Toeplitz matrix and B
* in the above format (related to A)
*
* given: A, B
* retrn: x (of Ax = B)
*
*-----------------------------------------------------------------------------
*/
MATRIX mat_lsolve_durbin( A, B )
MATRIX A, B;
{
MATRIX R, X;
int i, n;
n = MatRow(A);
R = mat_creat(n+1, 1, UNDEFINED);
for (i=0; i<n; i++)
{
R[i][0] = A[i][0];
}
R[n][0] = B[n-1][0];
X = mat_durbin( R );
mat_free( R );
return (X);
}


46
matrix/materr.c Normal file
View File

@ -0,0 +1,46 @@
/*
*-----------------------------------------------------------------------------
* file: materr.c
* desc: matrix error handler
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#ifdef __TURBOC__
#include <alloc.h>
#else
#include <malloc.h>
#endif
#include "matrix.h"
MATRIX mat_error( errno )
int errno;
{
switch( errno )
{
case MAT_MALLOC:
fprintf(stderr, "mat: malloc error\n" );
break;
case MAT_FNOTOPEN:
fprintf(stderr, "mat: fileopen error\n" );
break;
case MAT_FNOTGETMAT:
fprintf(stderr, "fgetmat: matrix read error\n");
break;
}
return (NULL);
}


77
matrix/matinv.c Normal file
View File

@ -0,0 +1,77 @@
/*
*-----------------------------------------------------------------------------
* file: matinv.c
* desc: matrix inversion
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi: 14 may 92 v0.2
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include <math.h>
#ifdef __TURBOC__
#include <alloc.h>
#else
#include <malloc.h>
#endif
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_inv
* desct: find inverse of a matrix
* given: a = square matrix a
* retrn: square matrix Inverse(A)
* NULL = fails, singular matrix, or malloc() fails
*-----------------------------------------------------------------------------
*/
MATRIX mat_inv( a )
MATRIX a;
{
MATRIX A, B, C, P;
int i, j, n;
double temp;
n = MatCol(a);
A = mat_copy(a);
B = mat_creat( n, 1, UNDEFINED );
C = mat_creat( n, n, UNDEFINED );
P = mat_creat( n, 1, UNDEFINED );
/*
* - LU-decomposition -
* also check for singular matrix
*/
if (mat_lu(A, P) == -1)
{
mat_free(A);
mat_free(B);
mat_free(C);
mat_free(P);
return (NULL);
}
for (i=0; i<n; i++)
{
mat_fill(B, ZERO_MATRIX);
B[i][0] = 1.0;
mat_backsubs1( A, B, C, P, i );
}
mat_free(A);
mat_free(B);
mat_free(P);
return (C);
}


60
matrix/matmul.c Normal file
View File

@ -0,0 +1,60 @@
/*
*-----------------------------------------------------------------------------
* file: matmul.c
* desc: matrix multiplication
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_mul
* desct: multiplication of two matrice
* given: A, B = compatible matrice to be multiplied
* retrn: NULL if malloc() fails
* else allocated matrix of A * B
* comen:
*-----------------------------------------------------------------------------
*/
MATRIX mat_mul( A, B )
MATRIX A, B;
{
int i, j, k;
MATRIX C;
if ((C = mat_creat( MatRow(A), MatCol(B), UNDEFINED )) == NULL)
return (NULL);
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(B); j++)
for (k=0, C[i][j]=0.0; k<MatCol(A); k++)
{
C[i][j] += A[i][k] * B[k][j];
}
return (C);
}
double mat_diagmul( A )
MATRIX A;
{
int i;
double result = 1.0;
for (i=0; i<MatRow(A); i++)
{
result *= A[i][i];
}
return (result);
}


136
matrix/matrix.h Normal file
View File

@ -0,0 +1,136 @@
/*
*-----------------------------------------------------------------------------
* file: matrix.h
* desc: matrix mathematics header file
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1b
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#ifndef SHUMATRIX
#define SHUMATRIX
/*
*-----------------------------------------------------------------------------
* internal matrix structure
*-----------------------------------------------------------------------------
*/
typedef struct {
int row;
int col;
} MATHEAD;
typedef struct {
MATHEAD head;
/*
* only the starting address of the following will be
* returned to the C programmer, like malloc() concept
*/
double *matrix;
} MATBODY;
typedef double **MATRIX;
#define Mathead(a) ((MATHEAD *)((MATHEAD *)(a) - 1))
#define MatRow(a) (Mathead(a)->row)
#define MatCol(a) (Mathead(a)->col)
/*
*----------------------------------------------------------------------------
* mat_errors definitions
*----------------------------------------------------------------------------
*/
#define MAT_MALLOC 1
#define MAT_FNOTOPEN 2
#define MAT_FNOTGETMAT 3
/*
*----------------------------------------------------------------------------
* matrice types
*----------------------------------------------------------------------------
*/
#define UNDEFINED -1
#define ZERO_MATRIX 0
#define UNIT_MATRIX 1
/* prototypes of matrix package */
#ifndef NOPROTO
MATRIX mat_error (int);
MATRIX _mat_creat (int, int);
MATRIX mat_creat (int, int, int);
MATRIX mat_fill (MATRIX, int);
int mat_free (MATRIX);
MATRIX mat_copy (MATRIX);
MATRIX mat_colcopy1 (MATRIX, MATRIX, int, int);
int fgetmat (MATRIX, FILE *);
MATRIX mat_dump (MATRIX);
MATRIX mat_dumpf (MATRIX, char *);
MATRIX mat_fdump (MATRIX, FILE *);
MATRIX mat_fdumpf (MATRIX, char *, FILE *);
MATRIX mat_add (MATRIX, MATRIX);
MATRIX mat_sub (MATRIX, MATRIX);
MATRIX mat_mul (MATRIX, MATRIX);
double mat_diagmul (MATRIX);
MATRIX mat_tran (MATRIX);
MATRIX mat_inv (MATRIX);
MATRIX mat_SymToeplz (MATRIX);
int mat_lu (MATRIX, MATRIX);
MATRIX mat_backsubs1 (MATRIX, MATRIX, MATRIX, MATRIX, int);
MATRIX mat_lsolve (MATRIX, MATRIX);
MATRIX mat_submat (MATRIX, int, int);
double mat_cofact (MATRIX, int, int);
double mat_det (MATRIX);
double mat_minor (MATRIX, int, int);
MATRIX mat_durbin (MATRIX);
MATRIX mat_lsolve_durbin(MATRIX, MATRIX);
#else
MATRIX mat_error ();
MATRIX _mat_creat ();
MATRIX mat_creat ();
MATRIX mat_fill ();
int mat_free ();
MATRIX mat_copy ();
MATRIX mat_colcopy1 ();
int fgetmat ();
MATRIX mat_dumpf ();
MATRIX mat_dump ();
MATRIX mat_fdump ();
MATRIX mat_fdumpf ();
MATRIX mat_add ();
MATRIX mat_sub ();
MATRIX mat_mul ();
double mat_diagmul ();
MATRIX mat_tran ();
MATRIX mat_inv ();
MATRIX mat_SymToeplz ();
int mat_lu ();
MATRIX mat_backsubs1 ();
MATRIX mat_lsolve ();
MATRIX mat_submat ();
double mat_cofact ();
double mat_det ();
double mat_minor ();
MATRIX mat_durbin ();
MATRIX mat_lsolve_durbin();
#endif
#endif

182
matrix/matsolve.c Normal file
View File

@ -0,0 +1,182 @@
/*
*-----------------------------------------------------------------------------
* file: matsolve.c
* desc: solve linear equations
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi: 14 may 92 v0.2
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
* [2] Kendall E.Atkinson, "An Introduction to Numberical Analysis,"
* John Wiley & Sons, 1978.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include <math.h>
#ifdef __TURBOC__
#include <alloc.h>
#else
#include <malloc.h>
#endif
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_lu
* desct: in-place LU decomposition with partial pivoting
* given: !! A = square matrix (n x n) !ATTENTION! see commen
* P = permutation vector (n x 1)
* retrn: number of permutation performed
* -1 means suspected singular matrix
* comen: A will be overwritten to be a LU-composite matrix
*
* note: the LU decomposed may NOT be equal to the LU of
* the orignal matrix a. But equal to the LU of the
* rows interchanged matrix.
*-----------------------------------------------------------------------------
*/
int mat_lu( A, P )
MATRIX A;
MATRIX P;
{
int i, j, k, n;
int maxi, tmp;
double c, c1;
int p;
n = MatCol(A);
for (p=0,i=0; i<n; i++)
{
P[i][0] = i;
}
for (k=0; k<n; k++)
{
/*
* --- partial pivoting ---
*/
for (i=k, maxi=k, c=0.0; i<n; i++)
{
c1 = fabs( A[(int)P[i][0]][k] );
if (c1 > c)
{
c = c1;
maxi = i;
}
}
/*
* row exchange, update permutation vector
*/
if (k != maxi)
{
p++;
tmp = P[k][0];
P[k][0] = P[maxi][0];
P[maxi][0] = tmp;
}
/*
* suspected singular matrix
*/
if ( A[(int)P[k][0]][k] == 0.0 )
return (-1);
for (i=k+1; i<n; i++)
{
/*
* --- calculate m(i,j) ---
*/
A[(int)P[i][0]][k] = A[(int)P[i][0]][k] / A[(int)P[k][0]][k];
/*
* --- elimination ---
*/
for (j=k+1; j<n; j++)
{
A[(int)P[i][0]][j] -= A[(int)P[i][0]][k] * A[(int)P[k][0]][j];
}
}
}
return (p);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_backsubs1
* desct: back substitution
* given: A = square matrix A (LU composite)
* !! B = column matrix B (attention!, see comen)
* !! X = place to put the result of X
* P = Permutation vector (after calling mat_lu)
* xcol = column of x to put the result
* retrn: column matrix X (of AX = B)
* comen: B will be overwritten
*-----------------------------------------------------------------------------
*/
MATRIX mat_backsubs1( A, B, X, P, xcol )
MATRIX A, B, X, P;
int xcol;
{
int i, j, k, n;
double sum;
n = MatCol(A);
for (k=0; k<n; k++)
{
for (i=k+1; i<n; i++)
B[(int)P[i][0]][0] -= A[(int)P[i][0]][k] * B[(int)P[k][0]][0];
}
X[n-1][xcol] = B[(int)P[n-1][0]][0] / A[(int)P[n-1][0]][n-1];
for (k=n-2; k>=0; k--)
{
sum = 0.0;
for (j=k+1; j<n; j++)
{
sum += A[(int)P[k][0]][j] * X[j][xcol];
}
X[k][xcol] = (B[(int)P[k][0]][0] - sum) / A[(int)P[k][0]][k];
}
return (X);
}
/*
*-----------------------------------------------------------------------------
* funct: mat_lsolve
* desct: solve linear equations
* given: a = square matrix A
* b = column matrix B
* retrn: column matrix X (of AX = B)
*-----------------------------------------------------------------------------
*/
MATRIX mat_lsolve( a, b )
MATRIX a, b;
{
MATRIX A, B, X, P;
int i, n;
double temp;
n = MatCol(a);
A = mat_copy(a);
B = mat_copy(b);
X = mat_creat(n, 1, ZERO_MATRIX);
P = mat_creat(n, 1, UNDEFINED);
mat_lu( A, P );
mat_backsubs1( A, B, X, P, 0 );
mat_free(A);
mat_free(B);
mat_free(P);
return (X);
}

42
matrix/matsub.c Normal file
View File

@ -0,0 +1,42 @@
/*
*-----------------------------------------------------------------------------
* file: matsub.c
* desc: matrix substraction
* by: ko shu pui, patrick
* date: 24 nov 91 v0.1
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_sub
* desct: subtraction of two matrice
* given: A, B = compatible matrice to be added
* retrn: NULL if malloc() fails
* else allocated matrix of A - B
* comen:
*-----------------------------------------------------------------------------
*/
MATRIX mat_sub( A, B )
MATRIX A, B;
{
int i, j;
MATRIX C;
if ((C = mat_creat( MatRow(A), MatCol(A), UNDEFINED )) == NULL)
return (NULL);
for (i=0; i<MatRow(A); i++)
for (j=0; j<MatCol(A); j++)
{
C[i][j] = A[i][j] - B[i][j];
}
return (C);
}

48
matrix/matsubx.c Normal file
View File

@ -0,0 +1,48 @@
/*
*-----------------------------------------------------------------------------
* file: matsubx.c
* desc: find submatrix
* by: ko shu pui, patrick
* date: 24 may 92 v0.4
* revi:
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_submat
* desct: return a submatrix S of A
* given: A = main matrix,
* i,j = row and column of A to be deleted to obtained S
* retrn: S
*-----------------------------------------------------------------------------
*/
MATRIX mat_submat( A, i, j )
MATRIX A;
int i,j;
{
int m, m1, p, p1;
MATRIX S;
S = mat_creat(MatRow(A)-1, MatCol(A)-1, UNDEFINED);
for (m=m1=0; m<MatRow(A); m++)
{
if (m==i) continue;
for (p=p1=0; p<MatCol(A); p++)
{
if (p==j) continue;
S[m1][p1] = A[m][p];
p1++;
}
m1++;
}
return (S);
}

43
matrix/mattoepz.c Normal file
View File

@ -0,0 +1,43 @@
/*
*-----------------------------------------------------------------------------
* file: mattoepz.c
* desc: matrix mathematics - toeplitz matrix
* by: ko shu pui, patrick
* date: 26 nov 91 v0.1
* revi: 14 may 92 v0.2
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_SymToeplz
* desct: create a n x n symmetric Toeplitz matrix from
* a n x 1 correlation matrix
* given: R = correlation matrix (n x 1)
* retrn: the symmetric Toeplitz matrix
*-----------------------------------------------------------------------------
*/
MATRIX mat_SymToeplz( R )
MATRIX R;
{
int i, j, n;
MATRIX T;
n = MatRow(R);
T = mat_creat(n, n, UNDEFINED);
for (i=0; i<n; i++)
for (j=0; j<n; j++)
{
T[i][j] = R[abs(i-j)][0];
}
return (T);
}

45
matrix/mattran.c Normal file
View File

@ -0,0 +1,45 @@
/*
*-----------------------------------------------------------------------------
* file: mattran.c
* desc: matrix mathematics
* by: ko shu pui, patrick
* date: v0.1 - 24 nov 91
* revi: v0.2 - 14 may 92
* ref:
* [1] Mary L.Boas, "Mathematical Methods in the Physical Sciene,"
* John Wiley & Sons, 2nd Ed., 1983. Chap 3.
*
*-----------------------------------------------------------------------------
*/
#include <stdio.h>
#include "matrix.h"
/*
*-----------------------------------------------------------------------------
* funct: mat_tran
* desct: transpose of a matrix
* given: A = matrix A to be transposed
* retrn: allocated matrix for A^t
* comen:
*-----------------------------------------------------------------------------
*/
MATRIX mat_tran( A )
MATRIX A;
{
int i, j;
MATRIX At;
if ((At = mat_creat( MatCol(A), MatRow(A), UNDEFINED )) == NULL)
return (NULL);
/*
* Transposing ...
*/
for (i=0; i<MatCol(A); i++)
for (j=0; j<MatRow(A); j++)
{
At[i][j] = A[j][i];
}
return (At);
}