An improved Gauss algorithm of the decision of rarefied systems of the linear equations has been considered in article.

author: Beskrovniy Konstantin
information: www.robosoft.info
published: 18.12.2006

Introduction

In many areas of a science and technics, at the decision of various problems, there is a necessity of solution of systems of the linear equations which sometimes are strongly rarefied. As usually such systems of the equations have the big order. To find their decision (in a little foreseeable terms) probably only with application of computer facilities. For the solve of this equations it is possible to use direct methods (of Gauss, of Cramer), and confidants methods (of Zeidel, of gradient release, etc.). Direct methods give the exact decision (the error can be only because of an overflow of type sizes at transformations for limits of accuracy of types for a used platform), but usually demand a lot of iterations. The approached methods usually work faster, but can demand special entry conditions or any preparation of the initial data, and the decision not always converges. In given article has review the realization of Gauss method modified for the decision of rarefied systems of the linear equations. The given method is not "mathematical", and more likely "algorithmic" and shows, how by simple receptions it is possible to speed up calculation in hundreds times.

Using conditions

The considered method approaches for anyone non confluent (having one solution) systems of the linear equations. The unique condition - a method is actual only for rarefied matrixes of big dimensions N> 1000 where not less than two thirds of coefficientes are equal 0. For usual matrixes the method will not give a prize neither in speed, nor in economy of memory.

Description of a method

The Gauss method consists in zeroing bottom (or top) half-matrix a[i,k], where row i=2..N, column k=1..i-1. Using a known rule, that the solution of system will not change, if to one of its equations to add another, increased on any coefficient, rows of a matrix serially, since the first, multiply on special coefficient and are summarized with all below lower lines. Coefficientes are selected so that each row is zeroing corresponding to it (i=k) column (only it under diagonal part). For more details about Gauss method click here.

The described algorithm is easily sold. The quantity of iterations is known beforehand and depends only on the order of a matrix. If him to program in classical style it will be easy to cope with small matrixes N~1000. However matrixes sometimes are about 100 000 and more. On calculation of such matrix by classical Gauss method there will leave years as the quantity of iterations will be huge. If the given matrix to store in a usual bidimentional array at use such as double it will occupy in memory 100 000 × 100 000 × 8 = 800 000 000 bytes = 800 Mb. And if to take into account that the size of a matrix grows proportionally to a square about a matrix, for a matrix about 400 000 it will be in 16 times more, i.e. 12.8 Gb!

Taking into account all mentioned facts it is possible to find out two directions of optimization of calculation:

• More compact storage of the information
• Reduction of quantity of iterations

For rarefied matrixes probably both variants:

More compact storage of the information - rarefied matrixes can be stored in memory not as a bidimentional array, and as some structure in which nonzero values of coefficientes are kept only, and also their coordinates in a matrix. Really, why store huge quantity of zero?. Reduction of quantity of iterations. If this structure to organize as special, that from it was possible to receive quickly the necessary information in the kind necessary for Gauss method, it is possible to reduce quantity of iterations up to minimal (in fact disappears necessity to touch zero coefficientes which in a matrix there can be 99 %).

For Gauss method would be conveniently fast to receive rows more to the right of the main diagonal (which will be multiplied and to be summarized with underlaying rows for zeroing bottom half-matrix), columns below main diagonal (actually that will be nulled) and elements of the main diagonal (values with which will be nulled columns). Therefore it is convenient to break a matrix on top half-matrix (everything, that is higher than the main diagonal), bottom half-matrix (everything, that is lower) and elements of the main diagonal. Thus it is necessary to transform the coordinate system:

• For the main diagonal it is one coordinate d which is equal to corresponding index of an element of the main diagonal
• For top and bottom half-matrix coordinates d, j, where j - displacement from the main diagonal to the right or down (as shown in figure)

Class diagram

Apparently from the diagram, communication of the basic components is carried out through interfaces, therefore the user can easily add the realizations of classes, not changing the general structure.

Realization of the method

All source codes are given in C#, and represent a full set of the classes necessary for functioning of a method. The code is written maximum simply for understanding of a principle of work of a method and consequently does not apply for the maximal productivity and is not optimum. About necessary completions it is told in the conclusion. The code is supplied with comments and does not require an additional explanation.

Interfaces of row and column:

public interface IRow
{
      // set the value in cell with index num
      void setValue(int num, double value);
	  
      // add value to existing cells value
      void addValue(int num, double value);

      // get value from cell with index num
      double getValue(int num);
	  
      // return all non-zero values of row/column:
      // indexes of cells - in indexes, values of cells in values
      void getValues(ref int[] indexes, ref double[] values);
}

Interface of matrix:

public interface IMatrix
{
    // set the value in cell with coordinates [row,col];
	// row – row number
	// col – column number
    void setValue(int row, int col, double value);

    // add value to existing value in cell with coordinates [row,col]
    void addValue(int row, int col, double value);

    // get value from cell with coordinates [row,col]
    double getValue(int row, int col);

    // get order of matrix
	int getN();

    // get all non-zero values from row with index d,
	// who located to right from main diagonal
    void getJRow(int d, ref int[] indexes, ref double[] values);
 
    // get all non-zero values from col with index d
	// who located to bottom from main diagonal
	void getJCol(int d, ref int[] indexes, ref double[] values);
}

Abstract class for a matrix. Includes fields and methods identical to different realizations of matrix class:

public abstract class Matrix : IMatrix
{
      // order of matrix
      private int N;
	  
      public Matrix(int N)
      {
            this.N = N;
      }
      // return order of matrix
      public int getN()
      {
            return N;
      }    

      // support IMatrix interface
      public abstract void setValue(int row, int col, double value);
      public abstract void addValue(int row, int col, double value);
      public abstract double getValue(int row, int col);
      public abstract void getJRow(int d, ref int[] indexes, ref double[] values);
      public abstract void getJCol(int d, ref int[] indexes, ref double[] values); 

      // service methods for conversing coordinates

      // return 0 for x0
      protected static int sigma(int x)
      {
            return (x>0)?x:0;
      }

      // converce row,col in half-matrix coordinates (d,j)
      // for upper half-matrix: d=0..(n-1), j=-(n-d-1)..-1
      // for bottom half-matrix:  d=0..(n-1), j=1..(n-d-1)
      protected static int getD(int row, int col)
      {
            return row-sigma(row-col);
      }

      protected static int getJ(int row, int col)
      {
            return col-row;
      }
 
      // opposite converce coordinates from (d,j) to (row,col)
      protected static int getRow(int d, int j)
      {
            return d + sigma(-j);
      }

      protected static int getCol(int d, int j)
      {
            return d + sigma(j);
      }
}

the given stage it is necessary to decide how to store values of cells of a matrix. For very strong rarefied matrixes when, in each row of a matrix of nonzero values it is not enough, the data are better for storing in usual arrays as search in a small array works quickly. If a matrix not strongly rarefied it is better to use the hash-table. Realizations of interface IRow for both cases are below resulted:

IRow interface realization with array use for storing

public class Array : IRow
{
      // array for indexes of non-zero cell
      private ArrayList indexes = new ArrayList();

      // array for values of non-zero cells
      private ArrayList values = new ArrayList();

      public void setValue(int num, double value)
      {
            int i = indexes.IndexOf(num); 
            // if zero in set, then remove this cell
            if (value==0)
            {
                  if (i!=-1)
                  {
                        indexes.RemoveAt(i);
                        values.RemoveAt(i);
                  }
                  return;
            }
                 
            // if non zero in set, then rewrite or add cell
            if (i!=-1)
            {
                  values[i] = value;
            }
            else
            {
                  indexes.Add(num);
                  values.Add(value);
            }
      }
 
     public void addValue(int num, double value)
      {
            int i = indexes.IndexOf(num);
            if (i!=-1)
            {
                  values[i] = (double)values[i] + value;
            }
            else
            {
                  indexes.Add(num);
                  values.Add(value);
            }
      }

      public double getValue(int num)
      {
            int i = indexes.IndexOf(num);
            if (i!=-1) return (double)values[i];
            return 0;
      }
 
      public void getValues(ref int[] indexes, ref double[] values)
      {
         indexes = (int[])this.indexes.ToArray(System.Type.GetType("System.Int32"));
         values = (double[])this.values.ToArray(System.Type.GetType("System.Double"));
      }
}

IRow interface realization with hash table use for storing:

public class  Hash : IRow
{
      // hash table in format: index_non_zero_cell = cell_value
      private Hashtable hash = new Hashtable();

      public void setValue(int num, double value)
      {
            // if zero in set, remove this cell
            if (value==0)
            {
                  if (hash.ContainsKey(num)) hash.Remove(num);
                  return;
            }
 
            // if non zero in set, then rewrite or add cell
            hash[num] = value;
      }
 
      public void addValue(int num, double value)
      {
            if (hash.ContainsKey(num)) hash[num] = (double)hash[num] + value;
              else hash[num] = value;
      }

      public double getValue(int num)
      {
            if (hash.ContainsKey(num)) return (double)hash[num];
            return 0;
      }

      public void getValues(ref int[] indexes, ref double[] values)
      {
            System.Collections.ICollection keys = hash.Keys;
            int count = keys.Count;
            indexes = new int[count];
            values = new double[count];
            int i=0;
            for (IEnumerator it = keys.GetEnumerator(); it.MoveNext();)
            {
                  indexes[i] = (int)it.Current;
                  values[i++] = (double)hash[it.Current];
            }
      }
}

Full IMatrix interface realization:

public class HashMatrix : Matrix
{
      // the store type: array or hash table
      public const byte ARRAY_TYPE = 0;
      public const byte HASH_TYPE = 1;
           
      // jp – cells of upper half-matrix
      // jm – cells of bottom half-matrix
      // dd – cells of main diagonal
      private IRow[] jp,jm;
      private double[] dd;
           
      // constructor
      // type – store type: ARRAY_TYPE/HASH_TYPE
      // N – order of matrix
      public HashMatrix(byte type, int N) : base(N)
      {
            dd = new double[N];
            if (type==HASH_TYPE)
            {
                  jp = new Hash[N];
                  jm = new Hash[N];
                  for (int i=0; inew Hash();
                             jm[i] = new Hash();
                        }
            }
            else
            {
                  jp = new Array[N];
                  jm = new Array[N];
                  for (int i=0; inew Array();
                        jm[i] = new Array();
                  }
            }
      }
 
      public override void setValue(int row, int col, double value)
      {
            if (row==col)
            {
                  dd[row] = value;
                  return;
            }
            int d = getD(row,col);
            int j = getJ(row,col);
            if (j>0) jp[d].setValue(j,value);
                  else jm[d].setValue(-j,value);
      }
 
      public override void addValue(int row, int col, double value)
      {
            if (value==0) return;
            if (row==col)
            {
                  dd[row] += value;
                  return;
            }
            int  d = getD(row,col);
            int  j = getJ(row,col);
            if  (j>0) jp[d].addValue(j, value);
                  else  jm[d].addValue(-j,value);
      }
 
      public override double  getValue(int  row, int  col)
      {
            if  (row==col) return  dd[row];
            int  d = getD(row,col);
            int  j = getJ(row,col);
            return  (j>0)?jp[d].getValue(j):jm[d].getValue(-j);
      }
 
      public override void  getJRow(int  d, ref int [] indexes, ref double [] values)
      {
            jp[d].getValues(ref  indexes, ref  values);
            for  (int  i=0; ipublic override void  getJCol(int  d, ref int [] indexes, ref double [] values)
      {
            jm[d].getValues(ref indexes, ref values);
            for  (int  i=0; i
And, at last, class for Gauss method implementation:
public class  Gauss
{
      // matrix 
      private IMatrix matrix; 

      // constructor, matrix – matrix with setted cells 
      public  Gauss(IMatrix matrix)
      {
            this .matrix = matrix;
      }
 
      // main void, return solution, accept the vector of free members 
	  public double [] calculate(double [] B)
      {
            // zeroing bottom half-matrix, take every each rows from up to down 
            // and summing with all rows who located bottom 
            for  (int  row=0; row<(matrix.getN()-1); row++)
            {
                  // get all non-zero values of zeroing column 
                  int [] colIndexes = new int [0];
                  double [] colValues = new double [0];
                  matrix.getJCol(row,ref  colIndexes,ref  colValues);
 
                  // get indexes and values of row cells who located to right from 
                  // main diagonal (left cells already is zero) 
                  int [] rowIndexes = new int [0];
                  double [] rowValues = new double [0];
                  matrix.getJRow(row, ref  rowIndexes, ref  rowValues);

                  // get the value of main diagonal cell 
                  double  dd = matrix.getValue(row,row);
                  for  (int  i=0; i// solve zeroing coefficient 
                     double  k = colValues[i]/dd;
 
                     // zeroing column cell, 
                     matrix.setValue(colIndexes[i],row,0);

                     // summing rows 
                     for  (int  ii=0; ii// summing free members 
                     B[colIndexes[i]] -= B[row]*k;
                  }
            }
 
            // create vector of solution 
            double [] x = new double [matrix.getN()];

            // use reverse motion find solution 
            for  (int  row = (matrix.getN()-1); row>=0; row--)
            {
               double  e = 0;
               int [] indexes = new int [0];
               double [] values = new double [0];
               matrix.getJRow(row,ref  indexes,ref  values);
               for  (int  i=0; ireturn  x;
      }
}

Conclusion
In article has review the method of the solution of systems of the linear equations on Gauss algorithm , and the code realizing this algorithm for C#. Though the given code also is full, for use in the user application it requires the following completion:

It is necessary to make processing of exclusive situations who can throws exception, for example such as division on zero or an an index out of bounds for arrays.
To realize settlement algorithm in the separate thread. If to cause a method calculate from the basic UI a thread, it will stop reception of messages of an environment (« will hang up the program ») till the moment of end of calculation.
If it is necessary to calculate some systems of the linear equations which differ only the right part (free members) naturally there is no necessity each time to zeroing bottom half-matrix. For this purpose it is necessary to realize class Gauss so that it nulling a matrix has store a history of changes. After that for calculation of another solution it will be necessary to transform only a vector of the free members transmitted to a method calculate according to a history. The history can be stored in a separate matrix, or in bottom half-matrix of basic matrix, who is not used at reverse motion
If the given algorithm is planned to realize on a platform, where direct work with pointer is possible (Delphi, C ++) it is better to use whenever possible them, it will allow to carry out smaller quantity of transformations in some places.
For some matrixes is necessary of realized of a factorization algorithm for exception of hit of zero from the main diagonal cells.
In implementation of interface IRow, at creation of instanses of arrays and hashtables, in an obvious kind property capacity was not set. In practice optimum to use speed and to keep memory, it is necessary to pick up compromise value of capacity. So with increase capacity also increase speed and the memory of use, with reduction - on the contrary increases. Optimum value capacity depends on structure of a matrix and does not exceed a maximum quantity of nonzero elements in row or column of a matrix.

The given algorithm (taking into account the specified remarks), was used in the RoboCAD project for calculation of systems of the linear equations of a method of final elements in static calculations of designs, and also in dynamic calculations at integration of systems of the differential equations of fluctuations which was reduced to the repeated decision of systems of the linear equations.
The literature and references

Andrew Troelsen E. C# and a .NET platform. Library of the programmer - Spb.: Peter, 2004.
Technical documentation of MSDN.