15,356,336 members
Articles / Programming Languages / C
Article
Posted 21 Oct 2020

19.1K views
9 bookmarked

Matrix Library in C

Rate me:
A matrix library in the C language that is useful for primary platforms and large data
This article develops a matrix algebra library in the C language. The library includes both basic and advanced functions that are defined in this regard. The library is properly coded to be fast and efficient by employing appropriate mathematical algorithms. Furthermore, there is no memory leak in the implemented functions that will be very critical for large matrices. The library depends on basic libraries that result in its ability to be simply exploited in the basic platforms.

Introduction

Matrix algebra is frequently encountered in a wide group of engineering and industrial applications. Generally, matrix algebra plays a key role in various engineering fields as a primary and basic tool. Hence, it is necessary to develop an efficient matrix library that contains the basic and advanced matrix functions such as `sum`, multiply, and inverse.

The matrix library should be fast with an efficient computational burden. Since this library is primarily used in engineering applications, it is very necessary to devote an efficient computational burden. This purpose can be followed by utilizing the best well-derived mathematical algorithms that are proposed for the matrix functions. There is a wide range of algorithms for each matrix function for solving these matrix functions. This article exploits the best algorithm that is developed in the literature to satisfy the computational burden constraints.

The developed matrix library depends on some basic and elementary libraries including Math.h. Hence, the library has a major potential to be used and implemented in the basic platforms involving micro-controllers. The library can be easily used which is explained later in this article.

The matrix library briefly contains a group of functions and structures. The functions are divided into two parts, namely elementary and advanced functions based on their background mathematical concepts. The elementary functions are a set of basic functions such as matrix multiplication and summation. The advanced functions are more complicated than the elementary ones which can be used in the advanced applications. All these functions are fully described in this article.

Background

Matrix algebra includes a set of matrix relations and functions that are widely used in engineering applications. This algebra is individually investigated in a lot of papers and books that basically describes its related topics and concepts. Some of these basic references are given below:

The above links provide some basic books about matrix algebra that can be reviewed to totally derive the results of this article. However, some of the implemented functions in this article are known and popular in various engineering fields.

A matrix is represented by a C `struct` in this article as shown here:

C++
```struct Mat{
double* entries;
int row;
int col;
};```

The above structure easily shows a matrix based on the idea of hashmap. This structure is repeatedly used as the input and output arguments of the implemented functions in the presented matrix library. As can be seen, the structure contains three fields that exactly represents the matrix.

Some developed functions have more than one output argument in the developed library. In this situation, another structure is exploited to return all final matrices which is a matrix list with the following definition:

C++
```struct MatList{
Mat* mat;
MatList* next;
};```

The implemented functions can be easily used by the above-mentioned structures and this point is totally explained in the next section.

Matrix Library

This article provides a matrix library that includes the matrix related functions that are categorized as elementary and advanced ones. First, the elementary functions are presented in the next subsection.

Elementary Functions

Some of these functions are provided to initial and define some specific matrices as given in the following code blocks. Functions `newmat()` and `freemat()` initializes and frees a matrix.

C++
```Mat* newmat(int r,int c,double d){
Mat* M=(Mat*)malloc(sizeof(Mat));
M->row=r;M->col=c;
M->entries=(double*)malloc(sizeof(double)*r*c);
int k=0;
for(int i=1;i<=M->row;i++){
for(int j=1;j<=M->col;j++){
M->entries[k++]=d;
}
}
return M;
}
void freemat(Mat* A){
free(A->entries);
free(A);
}```

There exist functions to initialize specific matrices that are presented as follows:

C++
```Mat* eye(int n){
Mat* I=newmat(n,n,0);
for(int i=1;i<=n;i++){
I->entries[(i-1)*n+i-1]=1;
}
return I;
}
Mat* zeros(int r,int c){
Mat* Z=newmat(r,c,0);
return Z;
}
Mat* ones(int r,int c){
Mat* O=newmat(r,c,1);
return O;
}```

Also, it is possible to initialize a matrix with random numbers with `randm()` as follows:

C++
```Mat* randm(int r,int c,double l,double u){
Mat* R=newmat(r,c,1);
int k=0;
for(int i=1;i<=r;i++){
for(int j=1;j<=c;j++){
double r=((double)rand())/((double)RAND_MAX);
R->entries[k++]=l+(u-l)*r;
}
}
return R;
}```

Above, parameters `r`and `c` declare the matrix that includes random numbers between real numbers `l` and `u`. This function is evaluated to initialize a 3 by 3 matrix whose entries lie within `-2` and `+2`. The result is presented in Figure 1 in the following:

Figure 1. The initialized random matrix.

The `set()` and `get()` functions are presented that will be used to get and set entry of the matrix located at a specific position:

C++
```double get(Mat* M,int r,int c){
double d=M->entries[(r-1)*M->col+c-1];
return d;
}
void set(Mat* M,int r,int c,double d){
M->entries[(r-1)*M->col+c-1]=d;
}```

Functions `sum()` and `minus()` compute the summation and minus of two pre-defined matrices:

C++
```Mat* sum(Mat* A,Mat* B){
int r=A->row;
int c=A->col;
Mat* C=newmat(r,c,0);
int k=0;
for(int i=0;i<r;i++){ for="">entries[k]=A->entries[k]+B->entries[k];
k+=1;
}
}
return C;
}

Mat* minus(Mat* A,Mat* B){
int r=A->row;
int c=A->col;
Mat* C=newmat(r,c,0);
int k=0;
for(int i=0;i<r;i++){ for="">entries[k]=A->entries[k]-B->entries[k];
k+=1;
}
}
return C;
}```

Functions `innermultiply()`, `scalermultiply()`, and `multiply()` compute the scaler and matrix multiplications:

C++
```double innermultiply(Mat* a,Mat* b){
double d=0;
int n=a->row;
if(a->col>n){
n=a->col;
}
for(int i=1;i<=n;i++){
d+=a->entries[i-1]*b->entries[i-1];
}
return d;
}

Mat* scalermultiply(Mat* M,double c){
Mat* B=newmat(M->row,M->col,0);
int k=0;
for(int i=0;i<m->row;i++){
for(int j=0;j<m->col;j++){
B->entries[k]=M->entries[k]*c;
k+=1;
}
}
return B;
}

Mat* multiply(Mat* A,Mat* B){
int r1=A->row;
int r2=B->row;
int c1=A->col;
int c2=B->col;
if (r1==1&&c1==1){
Mat* C=scalermultiply(B,A->entries[0]);
return C;
}else if (r2==1&&c2==1){
Mat* C=scalermultiply(A,B->entries[0]);
return C;
}
Mat* C=newmat(r1,c2,0);
for(int i=1;i<=r1;i++){
for(int j=1;j<=c2;j++){
double de=0;
for(int k=1;k<=r2;k++){
de+=A->entries[(i-1)*A->col+k-1]*B->entries[(k-1)*B->col+j-1];
}
C->entries[(i-1)*C->col+j-1]=de;
}
}
return C;
}```

The other useful elementary functions are used to manipulate some parts of a matrix. For example, `removerow()`and `removecolumn() `>are used for removing one specific row or column of the matrix, respectively.

C++
```Mat* removerow(Mat* A,int r){
Mat* B=newmat(A->row-1,A->col,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
if(i!=r){
B->entries[k]=A->entries[(i-1)*A->col+j-1];
k+=1;
}
}
}
return B;
}

Mat* removecol(Mat* A,int c){
Mat* B=newmat(A->row,A->col-1,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
if(j!=c){
B->entries[k]=A->entries[(i-1)*A->col+j-1];
k+=1;
}
}
}
return B;
}

Mat* submat(Mat* A,int r1,int r2,int c1,int c2){
Mat* B=newmat(r2-r1+1,c2-c1+1,0);
int k=0;
for(int i=r1;i<=r2;i++){
for(int j=c1;j<=c2;j++){
B->entries[k++]=A->entries[(i-1)*A->col+j-1];
}
}
return B;
}

Mat* copyvalue(Mat* A){
Mat* B=newmat(A->row,A->col,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
B->entries[k]=A->entries[k];
k++;
}
}
return B;
}```

The `transpose`, `trace`, determinant, `adjoint` of a matrix are obtained as given in the following code block:

C++
```Mat* transpose(Mat* A){
Mat* B=newmat(A->col,A->row,0);
int k=0;
for(int i=1;i<=A->col;i++){
for(int j=1;j<=A->row;j++){
B->entries[k]=A->entries[(j-1)*A->row+i-1];
k+=1;
}
}
return B;
}

double trace(Mat* A){
double d=0;
for(int i=1;i<=A->row;i++){
d+=A->entries[(i-1)*A->row+i-1];
}
return d;
}

double det(Mat* M){
int r=M->row;
int c=M->col;
if(r==1&&c==1){
double d=M->entries[0];
return d;
}
Mat* M1=removerow(M,1);
Mat* M2=newmat(M->row-1,M->col-1,0);
double d=0, si=+1;
for(int i=1;i<=M->col;i++){
double c=M->entries[i-1];
removecol2(M1,M2,i);
d+=si*det(M2)*c;
si=-si;
}
freemat(M1);
freemat(M2);
return d;
}

double trace(Mat* A){
double d=0;
for(int i=1;i<=A->row;i++){
d+=A->entries[(i-1)*A->row+i-1];
}
return d;
}

Mat* B=newmat(A->row,A->col,0);
Mat* A1=newmat(A->row-1,A->col,0);
Mat* A2=newmat(A->row-1,A->col-1,0);
for(int i=1;i<=A->row;i++){
removerow2(A,A1,i);
for(int j=1;j<=A->col;j++){
removecol2(A1,A2,j);
double si=pow(-1,(double)(i+j));
B->entries[(i-1)*B->col+j-1]=det(A2)*si;
}
}
Mat* C=transpose(B);
freemat(A1);
freemat(A2);
freemat(B);
return C;
}```

An example is provided for a better understanding of the above functions in the following:

Figure 2. A matrix and its related elementary results.

Finally, two types of inverse functions are developed which are ` inverse() ` and ` triinverse()`. These functions calculate the inverse matrix of a square matrix in general and triangular forms, respectively.

C++
```Mat* inverse(Mat* A){
double de=det(A);
Mat* C=scalermultiply(B,1/de);
freemat(B);
return C;
}

Mat* triinverse(Mat* A){
Mat* B=newmat(A->row,A->col,0);
for(int i=1;i<=B->row;i++){
for(int j=i;j<=B->col;j++){
if(i==j){
B->entries[(i-1)*B->col+j-1]=1/A->entries[(i-1)*A->col+j-1];
}else{
B->entries[(i-1)*B->col+j-1]=-A->entries[(i-1)*A->col+j-1]/A->
entries[(j-1)*A->col+j-1];
}
}
}
return B;
}   ```

This subsection presents advanced functions. First, a function `rowechelon()` is presented that obtains the row echelon form of a given matrix. The code is as follows:

C++
```Mat* rowechelon(Mat* A){
if(A->row==1){
for(int j=1;j<=A->col;j++){
if(A->entries[j-1]!=0){
Mat* B=scalermultiply(A,1/A->entries[j-1]);
return B;
}
}
Mat* B=newmat(1,A->col,0);
return B;
}
Mat* B=copyvalue(A);
int ind1=B->col+1;
int ind2=1;
for(int i=1;i<=B->row;i++){
for(int j=1;j<=B->col;j++){
if(B->entries[(i-1)*B->col+j-1]!=0&&j<ind1)
{ break="" ind1="j;" ind2="i;">1){
for(int j=1;j<=B->col;j++){
double temp=B->entries[j-1];
B->entries[j-1]=B->entries[(ind2-1)*B->col+j-1];
B->entries[(ind2-1)*B->col+j-1]=temp;
}
}
if(B->entries[0]!=0){
double coeff=B->entries[0];
for(int j=1;j<=B->col;j++){
B->entries[j-1]/=coeff;
}
for(int i=2;i<=B->row;i++){
coeff=B->entries[(i-1)*B->col];
for(int j=1;j<=B->col;j++){
B->entries[(i-1)*B->col+j-1]-=coeff*B->entries[j-1];
}
}
}else{
double coeff=0;
for(int j=1;j<=B->col;j++){
if(B->entries[j-1]!=0&&coeff==0){
coeff=B->entries[j-1];
B->entries[j-1]=1;
}else if (B->entries[j-1]!=0){
B->entries[j-1]/=coeff;
}
}
}
Mat* B1=removerow(B,1);
Mat* B2=removecol(B1,1);
Mat* Be=rowechelon(B2);
for(int i=1;i<=Be->row;i++){
for(int j=1;j<=Be->col;j++){
B->entries[i*B->col+j]=Be->entries[(i-1)*Be->col+j-1];
}
}
freemat(B1);
freemat(B2);
freemat(Be);
return B;
}```

The row echelon form is a useful function for computing some important properties of a matrix such as its `null` space. The row echelon form is basically obtained by a group of elementary Gaussian operations. An example is given for this function:

Figure 3. A matrix and its row echelon form.

As mentioned before, the row echelon form will be used the `null` space of a matrix that is obtained by the following function:

C++
```Mat* null(Mat *A){
Mat* RM=rowechelon(A);
int k=RM->row;
for(int i=RM->row;i>=1;i--){
bool flag=false;
for(int j=1;j<=RM->col;j++){
if(RM->entries[(i-1)*RM->col+j-1]!=0){
flag=true;
break;
}
}
if(flag){
k=i;
break;
}
}
Mat* RRM=submat(RM,1,k,1,RM->col);
freemat(RM);
int nn=RRM->col-RRM->row;
if(nn==0){
Mat* N=newmat(0,0,0);
return N;
}
Mat* R1=submat(RRM,1,RRM->row,1,RRM->row);
Mat* R2=submat(RRM,1,RRM->row,1+RRM->row,RRM->col);
freemat(RRM);
Mat* I=eye(nn);
Mat* T1=multiply(R2,I);
freemat(R2);
Mat* R3=scalermultiply(T1,-1);
freemat(T1);
Mat* T2=triinverse(R1);
freemat(R1);
Mat* X=multiply(T2,R3);
freemat(T2);
freemat(R3);
Mat* N=vconcat(X,I);
freemat(I);
freemat(X);
for(int j=1;j<=N->col;j++){
double de=0;
for(int i=1;i<=N->row;i++){
de+=N->entries[(i-1)*N->col+j-1]*N->entries[(i-1)*N->col+j-1];
}
de=sqrt(de);
for(int i=1;i<=N->row;i++){
N->entries[(i-1)*N->col+j-1]/=de;
}
}
return N;
}```

In the following, functions `vconcat()` and `hconcat()`are presented that concatenate two matrices in vertical or horizontal structures. The codes are given in the following:

C++
```Mat* hconcat(Mat* A,Mat* B){
Mat* C=newmat(A->row,A->col+B->col,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
C->entries[k]=A->entries[(i-1)*A->col+j-1];
k++;
}
for(int j=1;j<=B->col;j++){
C->entries[k]=B->entries[(i-1)*B->col+j-1];
k++;
}
}
return C;
}

Mat* vconcat(Mat* A,Mat* B){
Mat* C=newmat(A->row+B->row,A->col,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
C->entries[k]=A->entries[(i-1)*A->col+j-1];
k++;
}
}
for(int i=1;i<=B->row;i++){
for(int j=1;j<=B->col;j++){
C->entries[k]=B->entries[(i-1)*B->col+j-1];
k++;
}
}
return C;
}```

Matrix decomposition is another important topic that is implemented in this article. Two important matrix decompositions are implemented that are LU and OR decompositions. Their codes are separately given in the following code blocks. First, the LU decomposition is presented:

C++
```MatList* lu(Mat* A){
if(A->row==1){
MatList* ml=(MatList*)malloc(sizeof(MatList));
ml->mat=newmat(1,1,A->entries[0]);
ml->next=(MatList*)malloc(sizeof(MatList));
ml->next->mat=newmat(1,1,1);
return ml;
}
double a=A->entries[0];
double c=0;
if(a!=0){
c=1/a;
}
Mat* w=submat(A,1,1,2,A->col);
Mat* v=submat(A,2,A->row,1,1);
Mat* Ab=submat(A,2,A->row,2,A->col);
Mat* T1=multiply(v,w);
Mat* T2=scalermultiply(T1,-c);
Mat* T3=sum(Ab,T2);
MatList* mlb=lu(T3);
freemat(T1);
freemat(T2);
freemat(T3);
freemat(Ab);
Mat* L=newmat(A->row,A->col,0);
Mat* U=newmat(A->row,A->col,0);
int k=0;
for(int i=1;i<=A->row;i++){
for(int j=1;j<=A->col;j++){
if(i==1&&j==1){
L->entries[k]=1;
U->entries[k]=a;
k++;
}else if(i==1&&j>1){
U->entries[k]=w->entries[j-2];
k++;
}else if(i>1&&j==1){
L->entries[k]=c*v->entries[i-2];
k++;
}else{
L->entries[k]=mlb->mat->entries[(i-2)*mlb->mat->col+j-2];
U->entries[k]=mlb->next->mat->entries[(i-2)*mlb->next->mat->col+j-2];
k++;
}
}
}
MatList* ml=(MatList*)malloc(sizeof(MatList));
ml->mat=L;
ml->next=(MatList*)malloc(sizeof(MatList));;
ml->next->mat=U;
freemat(w);
freemat(v);
free(mlb);
return ml;
}  ```

Now, the QR decomposition is presented:

C++
```MatList* qr(Mat* A){
int r=A->row;
int c=A->col;
Mat* Q=newmat(r,r,0);
Mat* R=newmat(r,c,0);
Mat* ek=newmat(r,1,0);
Mat* uj=newmat(r,1,0);
Mat* aj=newmat(r,1,0);
for(int j=1;j<=r;j++){
submat2(A,aj,1,r,j,j);
for(int k=1;k<=r;k++){
uj->entries[k-1]=aj->entries[k-1];
}
for(int k=1;k<=j-1;k++){
submat2(Q,ek,1,r,k,k);
double proj=innermultiply(aj,ek);
for(int l=1;l<=ek->row;l++){
ek->entries[l-1]*=proj;

}
uj=minus(uj,ek);
}
double nuj=norm(uj);
for(int i=1;i<=r;i++){
Q->entries[(i-1)*r+j-1]=uj->entries[i-1]/nuj;
}
for(int j1=j;j1<=c;j1++){
R->entries[(j-1)*c+j1-1]=innermultiply(uj,submat(A,1,r,j1,j1))/nuj;
}
}
MatList* ml=(MatList*)malloc(sizeof(MatList));
ml->mat=Q;
ml->next=(MatList*)malloc(sizeof(MatList));;
ml->next->mat=R;
freemat(ek);
freemat(uj);
freemat(aj);
return ml;
}```

Finally, an example is considered for these matrix decomposition methods:

Figure 4. A random matrix and its decompositions.

In the above example, the main code is as follows:

C++
```int _tmain(int argc, _TCHAR* argv[])
{
Mat* A=randm(4,4,-5,5);
printf("A=");
showmat(A);
MatList* ml=lu(A);
Mat* L=ml->mat;
Mat* U=ml->next->mat;
MatList* ml1=qr(A);
Mat* Q=ml1->mat;
Mat* R=ml1->next->mat;
printf("LU decomposition: \n");
printf("L=");
showmat(L);
printf("U=");
showmat(U);
printf("QR decomposition: \n");
printf("Q=");
showmat(Q);
printf("R=");
showmat(R);
getchar();
return 0;
}```

Using the Code

A matrix library can be simply used by including its header file in your module as presented in the following:

C++
```//
#include "MatLib.h"
//```

Points of Interest

This article proposes a matrix library that has major properties. It is fast due to employing efficient mathematical algorithms. It does not use advanced libraries which enables us to employ it in simple platforms. Moreover, in this library, each function optimally collects its memory garbages.

History

• 21st October, 2020: Initial version
• 24st October, 2020: Second version

Share

 Software Developer Iran (Islamic Republic of)
No Biography provided

 First PrevNext
 there's a bug in transpose() julian1 202220-Jan-22 12:45 julian1 2022 20-Jan-22 12:45
 cant get triinverse to work correctly Member 117206811-Nov-20 6:22 Member 11720681 1-Nov-20 6:22
 Thank you Member 1495827528-Oct-20 13:48 Member 14958275 28-Oct-20 13:48
 Excellent Code Member 1172068126-Oct-20 10:12 Member 11720681 26-Oct-20 10:12
 Re: Excellent Code Roozbeh Abolpour27-Oct-20 8:18 Roozbeh Abolpour 27-Oct-20 8:18
 Re: Excellent Code Member 1172068128-Oct-20 7:05 Member 11720681 28-Oct-20 7:05
 Re: Excellent Code Roozbeh Abolpour30-Oct-20 0:26 Roozbeh Abolpour 30-Oct-20 0:26
 Re: Excellent Code Member 1172068130-Oct-20 4:26 Member 11720681 30-Oct-20 4:26
 broken link Southmountain24-Oct-20 9:05 Southmountain 24-Oct-20 9:05
 Re: broken link Roozbeh Abolpour24-Oct-20 9:13 Roozbeh Abolpour 24-Oct-20 9:13
 Re: broken link Southmountain26-Oct-20 5:48 Southmountain 26-Oct-20 5:48
 MatLib.rar works, but MatLib.zip not.. diligent hands rule....
 Error Detection JohnLewis194523-Oct-20 11:07 JohnLewis1945 23-Oct-20 11:07
 Re: Error Detection Roozbeh Abolpour24-Oct-20 3:07 Roozbeh Abolpour 24-Oct-20 3:07
 Possible memory and performance issues Daniele Rota Nodari22-Oct-20 22:12 Daniele Rota Nodari 22-Oct-20 22:12
 Re: Possible memory and performance issues Roozbeh Abolpour24-Oct-20 5:53 Roozbeh Abolpour 24-Oct-20 5:53
 Re: Possible memory and performance issues Southmountain24-Oct-20 8:53 Southmountain 24-Oct-20 8:53
 Another Suggestion Rick York21-Oct-20 5:11 Rick York 21-Oct-20 5:11
 Re: Another Suggestion Roozbeh Abolpour21-Oct-20 6:59 Roozbeh Abolpour 21-Oct-20 6:59
 Re: Another Suggestion Rick York21-Oct-20 7:54 Rick York 21-Oct-20 7:54
 Re: Another Suggestion Roozbeh Abolpour21-Oct-20 8:14 Roozbeh Abolpour 21-Oct-20 8:14
 Re: Another Suggestion Rick York21-Oct-20 14:48 Rick York 21-Oct-20 14:48
 Re: Another Suggestion Roozbeh Abolpour21-Oct-20 19:56 Roozbeh Abolpour 21-Oct-20 19:56
 Re: Another Suggestion Southmountain21-Oct-20 8:41 Southmountain 21-Oct-20 8:41
 Re: Another Suggestion Roozbeh Abolpour21-Oct-20 8:51 Roozbeh Abolpour 21-Oct-20 8:51
 Re: Another Suggestion Armando A Bouza24-Oct-20 3:59 Armando A Bouza 24-Oct-20 3:59
 Last Visit: 31-Dec-99 18:00     Last Update: 4-Jul-22 19:42 Refresh 12 Next ᐅ