-----------------------------------
matrices.h
-----------------------------------

"matrices.h" is an header file for C floating point double precision matrix computation. It implements the main operations on matrices: multiplication (row by column), sum of matrices, multiplication by a scalar, calculation of the transposed matrix, calculation of the inverse matrix and the determinant of a square matrix, solving linear systems.

It includes functions for dynamic allocation of matrices and for input/output of matrices from/to file and from/to keyboard/screen.
It requires the header files <stdio.h> and <stdlib.h> .

All the source code on this page is released under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version, and is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.

Table of contents

Description of the functions and their use

Each matrix is represented by a structure comprising a double pointer to double precision floating point data (**el) for the elements and two integers for the orders, one for the number of rows (ordi) and one for the number of columns (ordj).  Elements are stored by lines (arrays of double precision floating point cells), pointed to by a vector of addresses.

typedef struct matrix{
 double **el;
 int ordi;
 int ordj;}MAT;

You will declare matrices as "MAT" objects in your programs:
MAT a;
a.ordi and a.ordj are respectively the number of rows and columns of the matrix "a".
The element (i,j) of the matrix "a" is a.el[i][j] .


- The "allocm" function allocates memory for a matrix "a" of orders "ordi" and "ordj":
allocm(&a, ordi, ordj);


- The "freem" function frees the memory allocated for the matrix "a:
freem(&a);


- The "msum" function neatly adds the elements of the matrices "a1" and "a2":
msum(&sum, a1, a2);
It operates on matrices of the same orders (rows and columns), including rectangular ones.
Automatically allocates memory for the "sum" matrix.


- The "mmulsc" function multiplies all elements of a matrix "a" by the scalar "lambda":
mmulsc(&prod, a, lambda);
It also works on rectangular matrices.
Automatically allocates memory for the "prod" product matrix.


- The "mmult" function executes the scalar product between the rows of the first matrix (a1) and the columns of the second one (a2):
mmult(&prod, a1, a2);
It also operates on rectangular matrices with the constraint that the number of columns of the first matrix must be equal to the number of rows of the second one.
Automatically allocates memory for the "prod" product matrix.


- The "trasnp" function transposes the matrix "a":
transp(&t, a);
It also works on rectangular matrices.
Automatically allocates memory for the transposed matrix "t".


- The "det" function recursively calculates the determinant "d" of the matrix "a" with the Laplace method (expansion of the matrix through the cofactors):
d=det(a);
It operates exclusively on square matrices.
Returns the value of the determinant as a double-precision floating point number.


- The "inv" function calculates the inverse of the matrix "a":
inv(&in, a);
It operates exclusively on square matrices.
Automatically allocates memory for the "in" inverse matrix.


- The "linsys" function calculates the solutions of a linear system with Cramer's rule given the matrix of the coefficients of the unknowns "a" and the column vector of the constant terms "v":
linsys(&sol, a, v);
It operates exclusively on square matrices.
Automatically allocates memory for the solution column vector "sol".
The column vector of the constant terms "v" must have a number of elements equal to the order of the square matrix of the coefficients.


- The writem / fwritem functions write an "a" matrix to standard output or file respectively:
writem(a);
fwritem(file_handler, a);


- The readm / freadm functions read an "a" matrix from standard input or from file respectively:
readm(&a);
freadm(file_handler, &a);
They automatically allocate memory for the "a" matrix.

Example of use (1): solving a linear system of equations

The "syslin.c" program uses the functions declared in the header file "matrices.h" to solve a linear system of equations given the matrix of the coefficients of the unknowns and the vector of constant terms.
Suppose you want to solve the following linear system of 6 equations:

3a +2b +7c +11d -8e +f = 8
5a -9b +2c +7d +e -2f = -1
4a +2b -5c +6d +13e +9f = 7
2a +6b +9c -3d -4e -f = 2
7a +7b -4c -5d -6e +12f = -5
3a +9b +17c +6d -11e +4f = -2

First of all you must create two text files with a text editor, one for the matrix of the coefficients, let's say "coeff.txt" and the other for the vector of constant terms, let's say "vect.txt":
The first row of each file will contain the number of row and cols (integers); the other rows will contain the elements in floating point format.

So we will have the following text file "coeff.txt" for the matrix of coefficients:

6  6
3.0  2.0  7.0  11.0  -8.0  1.0
5.0  -9.0  2.0  7.0  1.0  -2.0
4.0  2.0  -5.0  6.0  13.0  9.0
2.0  6.0  9.0  -3.0  -4.0  -1.0
7.0  7.0  -4.0  -5.0  -6.0  12.0
3.0  9.0  17.0  6.0  -11.0  4.0

and the following text file "vect.txt" for the vector of constant terms:

6  1
8.0
-1.0
7.0
2.0
-5.0
-2.0

At this point, to solve the linear system you must invoke the program "syslin" as follows; let's say that your output file is "sol.txt":
syslin coeff.txt vect.txt sol.txt
The file "sol.txt", created by the program, contains the vector of solutions:

6  1
1.319879
1.975108
-1.099400
1.129452
0.320457
-2.074376

The solution of the system is therefore:
a = 1.319879
b = 1.975108
c = -1.099400
d = 1.129452
e = 0.320457
f = -2.074376

Compiling and executing syslin.c

To download "syslin.c", the header file "matrices.h" needed to compile it, the two example files "coeff.txt" and "vect.txt" and the output file "sol.txt" as produced with the given input go to the download section of this web page.

Example of use (2): inversion of a square matrix

The invers.c program inverts a square matrix and prints the inverse matrix on the screen.
Suppose you want to invert the following 4x4 floating point square matrix:

3  2  7  11
5  -9  22  7
4  12  -5  6
1  -6  9  -3

First of all create a text file for the square matrix to be inverted with a text editor, let's say "matrix.txt"; the first two numbers will be the number of rows and cols (integers), the others will be the elements (floating point).

4  4
3.0  2.0  7.0  11.0
5.0  -9.0  22.0  7.0
4.0  12.0  -5.0  6.0
1.0  -6.0  9.0  -3.0

After creating the input file "matrix.txt" invoke the program "invers" as follows:
invers matrix.txt
The inverse matrix will be displayed on the screen.
The program displays an error message if you attempt to invert a non square matrix or if the square matrix you want to invert is non invertible (singular).

Compiling and executing invers.c

To download "invers.c", the header file "matrices.h" needed to compile it and the example input file "matrix.txt" go to the download section of this web page.

Downloads

Header file "matrices.h":  matrices.h

Linear system solver:  syslin.c

Matrix of the coefficients for the linear system solver:  coeff.txt

Vector of the constant terms for the linear system solver:  vect.txt

Vector of the solutions as produced by the linear system solver:  sol.txt

Square matrix inversion:  invers.c

Matrix to be inverted:  matrix.txt

Bug fixes

06/01/2022 fixed a bug in function allocm (bad name of the pointer to the rows , was "pmat->m", now corrected to "pmat->el")
06/01/2022 fixed a bug in function freem (bad name of the pointer to the rows , was "pmat->m", now corrected to "pmat->el")