DTBX: Digital Toolbox (C/C++ API library)

This software can be redistributed under GNU Lesser General Public License.
Every source code of this software can be obtained through GitHub
Copylight (c) 2019 Shigeo Kobayashi. All rights reserved.

Digital Toolbox (C/C++ API library) offers some useful functions(as of 2019) such as:

Matrix block triangular decomposition
Linear equation system solver
Sorting

This software has been tested on Windows-10(32-bit&64-bit) and Linux(32-bit CentOS-5 & 64-bit CentOS-7).
To build this software,include "dtbx.h" in your program and

Linux ... see makefile bundled.
Windows ... follow Visual studio instructions.

Matrix block triangular decomposition

Matrix block triangular decomposition is sometimes useful to solve linear equation system such as

    Ax=y

,where A is a square matrix,x and y are solution and constant vector respectively.
The following explanation is also applicable to decompose a large non-linear equation system to a sequence of smaller sized non-linear equation systems.
In this case,A can be a Jacobian matrix of non-linear eqation system.

It is sometimes possible to decompose A into a sequence of smaller sized sub-matrices as

    Ax=y  ===>  A_ix_i=y_i i=0,1,2...

where A_i and x_i are subset of A and x respectively. But y_i is not a subset of y that can be derived from y,A and x.
For example,the left matrix A(and x,y) shown bellow can be rearranged to the right one by exchanging some columns and rows.
(1 means the element of any non-zero value.)

     1 0 0 1 0 0 1 1  x₀ y₀         A₀ 1 1 0 0 0 0 0 0  x₄ y₃
     1 0 0 0 0 1 1 1  x₁ y₁         A₀ 1 1 0 0 0 0 0 0  x₃ y₅
     0 0 0 0 0 1 0 1  x₂ y₂         A₁ 0 0 1 1 1 0 0 0  x₇ y₄
     0 0 0 1 1 0 0 0  x₃ y₃         A₁ 0 0 1 1 1 0 0 0  x₂ y₆
  A  0 1 1 0 0 0 0 1  x₄ y₄   ==>   A₁ 0 0 1 1 1 0 0 0  x₁ y₇
     0 0 0 1 1 0 0 0  x₅ y₅         A₂ 0 0 1 0 0 1 0 0  x₅ y₂
     0 1 1 0 0 0 0 1  x₆ y₆         A₃ 0 1 1 0 0 0 1 1  x₆ y₀
     0 1 1 0 0 0 0 1  x₇ y₇         A₃ 0 0 1 0 0 1 1 1  x₀ y₁

Seeing the right 8x8 matrix,it is obvious that Ax=y can be divided into 4 subsets as A_ix_j=y_k,i = 0,1,2,3
Where the size of A₀ is 2x2,A₁ is 3x3,A₂ is 1x1,and A₃ is 2x2.
This result can be obtained by running the test program included in this release.
See the results of DtxDivideMatrix() in the test program(test.c).

Ignoring the actual form of A, y can be represented as y_i=f_i(x_j,x_k,x_l,...). in general form.
For instance y₀=f₀(x₀,x₃,x₆,x₇),in above example(see row₀ of the left matrix).
In graphically,edges(arrows) starting from nodes x₀,x₃,x₆,and x₇ end to y₀.

The steps to decompose A are:

Make 1 to 1 correspondings between y and x(pairing). This is the same as to set all diagonal elements of A be non-zero by exchanging some rows and columns.
Add reverse arrow from each y to paired x. This creates small loops(total number of the loops is the size of A).
Merge loops in the same order to larger loops. Resulting loops construct equation systems.

Y's and x's in the same loop construct 1 equation system solved.
The order of loops are the computation order of equation systems resulted.

int DtxDivideMatrix(int n,void A,int y[],int x[],int sizes[],int (IsNonZeroElement)(void *A,int I,int J,int n))

DtxDivideMatrix() divides the matrix given(=A) into a sequence of smaller sized matrices.
Because the element value of the given matrix is zero or non-zero is the only information needed,this function is not interested in the type of matrix(and it's element type),the user must provide the function IsNonZeroElement() which tells the specified element A[I,J] is zero or non-zero.
The resulting divided submatrices can be obtained by y[],x[],and sizes[].

  n      ... [IN]  The size of the (square)matrix divided.
  A      ... [IN]  (square)matrix divided or any information for the matrix,which is passed to the user 
                   supplied function IsNonZeroElement(). DtxDivideMatrix() never touch this A.
  y[]    ... [OUT] Row numbers of A are stored on exit.
  x[]    ... [OUT] Column numbers of A are stored on exit.
  sizes[]... [OUT] Sizes of divided sub-matrices.
                   If sizes[0] is 3,then rows from y[0] to y[2] and columns from x[0] to x[2]
                   form the first divided sub-matrix of A. 
  IsNonZeroElement(void *A,int I,int J,int n)
         ... [IN]  The user specified function which indicates the element A[I,J] is zero(returns 0)
                   or non-zero(return non-zero value). n is the size of A.
                   Typically,IsNonZeroElement() would be:

                   int IsNonZeroElement(void *A,int i,int j,int n)
                   {
	                if(((double*)A)[INDEX2(i,j,n)]!=0.0) return 1;
	                return 0;
                   }


  returns > 0;
         ... totoal number of sub-matrices when A is successfully divided.
  other  ... A is not divided (A is a singular matrix).

This DtxDivideMatrix() internally calls DtxSelectPair() and DtxDividePairedMatrix()
The usage of DtxDivideMatrix() can be found in the source codes of DtxLeqDivSolve().

int DtxSelectPair(int n,int lpair[],int rpair[],int num_r[],int *rv[])

DtxSelectPair() creates 1 to 1 correspondences(pairs) between left(y) and right(x) hand side variables.

   n             ... [IN]  total number of variables(total number of y or x).
   paired_x[n]   ... [OUT] 1 to 1 correspondences resulted. 
                           the left hand side variable(y) i is related to the right hand side variable paired_x[i].
   paired_y[n]   ... [OUT] 1 to 1 correspondences resulted.
                           the right hand side variable(x) i is related to the left hand side variable paired_y[i].
   num_r[n]      ... [IN]  num_r[i] is the total number of the right hand side variables of the variable i.
   rv[]          ... [IN]  rv[i] is the pointer to the array holding the right hand side variables of i. 
                           rv[i][j] is the jth right hand side variable of the left hand side variable i.

   returns 0     ... All varables are successfully paired each other.
           other ... Failed to make pairs.

After successfull completion,y[i] is paired with x[paired_x[i]], or x[i] is paired with y[paired_y[i]].
num_r[] and rv[] are data structure describing the functional relations between y and x.
For example, if y₀=f₀(x₀,x₃,x₆,x₇) then

  num_r[0] = 4 
  rv[0]    = {0,3,6,7}

The results of this function can be passed to the next DtxDividePairedMatrix().
But,it must be noted that the contents of rv[] must be renumbered so that y_i is paired with x_i
because DtxDividePairedMatrix() assumes y_i is paired with x_i.
The example codes for renumbering of rv[] can be found in DtxDivideMatrix().

int DtxDividePairedMatrix(int n,int order[],int block[],int num_r[],int *rv[])

DtxDividePairedMatrix() divides a matrix(or equation system) into a series of smaller sized matrices.
This function assumes that the left hand side valiable y_i and the right hand side variable x_i are paired each other.

 
   n        ... [IN]  the total number of variables(total number of y or x).
   order[n] ... [OUT] variables sorted to the order of matrices devided.
   block[n] ... [OUT] sizes of devided matrices.
   num_r[n] ... [IN]  the number of the right hand side variables. The nuber of the right hand side variables of the variable i is num_r[i].
   rv[]     ... [IN]  pointer to the array having the right hand side variables. The right hand size variavles of the variavle i are stored to rv[i].

   returns the number of total matrices divied (>0),otherwise negative error code.

Linear equation system solver

Linear equation system solver offers functions to solve linear equation systems such as:

  Ax = b

where A is a coefficient matrix of size n,x is a solution vector,and b is a constant vector.
The main functions are:

int DtxLeqSolve(double x[], double A[], double b[], int n)
Solves Ax=b directly(after LU decomposition)
int DtxLeqDivSolve(double x[], double A[], double b[], int n)
Solves Ax=b by decomposing into smaller sized equations by calling DtxDivideMatrix() and DtxLeqSolve().

Where A is a coeficient matrix of size nxn,x and b are solution and constant vectors of size n respectively.
DtxLeqSolve() and DtxLeqDivSolve() returns 0 after successive computation,or returns negative value when A is singular or any error occures.
Both functions allocate some working memeries and free them on exit which may be inefficient if they are called frequently (it is easy to change such memories be prepaired by caller,see source codes.).
Both functions destroy the contents of A,A should be saved before calling if necessary.

[Matrix indexing]

Because C does not support multi-dimensional array like FORTRAN or other languages.
To access A_i,j element of A,use INDEX2(i,j,n) macro defined in dtbx.h as
  #define INDEX2(i,j,M)   ((i)*(M)+j)
If FORTRAN compatible memory allocation is prefferable,redefine INDEX2(i,j,M) as
  #define INDEX2(i,j,M)   ((j)*(M)+i)
and re-compile everything.

Linear equation system solver offers 2 more fuinctions:

int DtxLuDecompose(double A[], double scales[], int index[], int n)

   decompose given matrix A to LU form. L is the lower triangular, and U is the upper triangular matrix.
   L and U are stored A.
   A         ... [IN/OUT], Matrix to be decomposed,and result of decomposition are stored on exit.
   scales[n] ... [OUT], work array for pivotting.
   index [n] ... [OUT], indices for pivotting,which is used later in DtxLuSolve()

   returns 0,  normally decomposed.
          -k(negative), failed to decompose when processing A[k-1,*].

void DtxLuSolve(double x[], double A[], double b[], int index[], int n)

   solves linear equation system Ax=b, A must be processed by DtxLuDecompose() before hand.
   index[n] must be the array also processed by DtxLuDecompose().

Linear equation system solver by indexing

DtxLeqDivSolve() internally allocates sub-matrices obtained by dividing A.
The following 3 functions use indexing instead of allocating sub-matrices.
Arrays I[] and J[] specify indexes of sub-matrices.
n ... size of sub-matrix.
N ... size of A and x,b.

int DtxLeqDivIndexSolve(double x[], double A[], double b[], int n)
int DtxLuIndexDecomp(double A[],int N,int I[], int J[], int n) ... DtxLeqDivIndexSolve() prepaires I[] and J[].
void DtxLuIndexSolve(double x[], double LU[], double b[],int N,int I[],int J[],int n)

Note:A and b are [IN/OUT],x is [OUT].

Sorting

Following 3 functions are available for sorting.

  heap   sort: void DtxHeapSort  (void *Target,int n,int (*compare)(void *Target,int i,int j),void (*swap)(void *Target,int i,int j) )
  bubble sort: void DtxBubbleSort(void *Target,int n,int (*compare)(void *Target,int i,int j),void (*swap)(void *Target,int i,int j) )
  quick  sort: void DtxQuickSort (void *Target,int n,int (*compare)(void *Target,int i,int j),void (*swap)(void *Target,int i,int j) )

void *Target is a pointer of the target to be sorted.
These 3 functions are not interrested in the target type being sorted,because comparison and exchange of 2 elements are done by the user.
The last 2 arguments,compare and swap, perform actual sort operations(compares elements and exchanges(swaps) elements of void *Target respectively.
void *Target can be anything of which element can be identified by it's index.
n is the element size of Target.

int (*compare)(void *Target,int i,int j) must return:
     positive value(>0) if Target[i]>Target[j]
     0                  if Target[i]==Target[j]
     negative value(<0) if Target[i]<Target[j]

void (*swap)(void *Target,int i,int j) must exchnage Target[i] and Target[j].

Shigeo Kobayashi 2019-11-30