need help improving matrix calculator

[ronenp88]

Hello everyone, I coded recently this matrix calculator
https://github.com/ronenp/MatrixCalculator
It works fine for small matrices, but for very large ones, 80 by 80 and so,
I think variables overflow, even with long long int, its because of the way i handle the calculation,
to get accurate results I don't deal with fractions until the very last step where I am about to print the solution, I multiply equations to a LCM then subtract or add depend on signs, but this method gets to big numbers fast, again for small matrices its fine, but for large complex ones, this algorithm just can't manage it, how should i improve or change the way this calculator work on those matrices? any suggestions please?

[ronenp88]

1	1	1	1	1	1	1
11	9	14	5	-1	-1	38
4	10	12	2	2	7	48
3	7	12	3	4	1	48
8	11	3	-3	4	11	49
7	13	12	6	3	1	41

this matrix for example if you want to check out what happens,
if you remove the comments on lines: 122,123 and on 187, 188
which will print to the console the process, you will see it get stuck when overflow occurs, this is pity
anyone knows how to solve this problem? maybe some workaround?

[lastchance]

@ronenp88,
what are you actually trying to do? "Matrix calculator" means WHAT?

If you leave the result of adding several unrelated fractions in rational form then your lowest common denominator will eventually overflow any intrinsic type (long long or otherwise). So it's a waste of time. Use doubles.

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <cassert>
using namespace std;

const double SMALL = 1.0E-30;
const double NEARZERO = 1.0E-10;
bool doPivot = true;

using vec    = vector<double>;
using matrix = vector<vec>;

void checkGauss();
void checkGaussInversion();
void print( const string &title, const matrix &A );
void print( const string &title, const vec &V );
matrix matmul( const matrix &A, const matrix &B );
vec matmul( const matrix &A, const vec &V );
matrix identity( int n );
bool GaussElimination( matrix A, vec B, vec &X );
bool GaussInversion( matrix A, matrix &B );

//===

int main()
{
   checkGauss();
   checkGaussInversion();
}

//===

void checkGauss()
{
   matrix A = { {  1 , 1 , 1 , 1 , 1 , 1  },
                {  11, 9 , 14, 5 , -1, -1 },
                {  4 , 10, 12, 2 , 2 , 7  },
                {  3 , 7 , 12, 3 , 4 , 1  }, 
                {  8 , 11, 3 , -3, 4 , 11 },  
                {  7 , 13, 12, 6 , 3 , 1  } };
   vec B = { 1, 38, 48, 48, 49, 41 };

   cout << "\n\n== Gaussian elimination to solve AX = B ===\n";
   vec X;
   if ( GaussElimination( A, B, X ) )
   {
      print( "X:" , X );
      print( "B:" , B );
      print( "AX:", matmul( A, X ) );
   }
   else
   {
      cout << "Unable to solve\n";
   }
}

//===

void checkGaussInversion()
{
   matrix A = { {  1 , 1 , 1 , 1 , 1 , 1  },
                {  11, 9 , 14, 5 , -1, -1 },
                {  4 , 10, 12, 2 , 2 , 7  },
                {  3 , 7 , 12, 3 , 4 , 1  }, 
                {  8 , 11, 3 , -3, 4 , 11 },  
                {  7 , 13, 12, 6 , 3 , 1  } };
   matrix B;

   cout << "\n\n== Calculate B = Ainverse ==\n";
   if ( GaussInversion( A, B ) )
   {
      print( "A:" , A );
      print( "B:" , B );
      print( "AB:", matmul( A, B ) );
   }
   else
   {
      cout << "Unable to invert\n";
   }
}

//===

void print( const string &title, const matrix &A )
{
   if ( title != "" ) cout << title << '\n';

   for ( auto &row : A )
   {
      for ( auto x : row )
      {
         if ( abs( x ) < NEARZERO ) x = 0.0;
         cout << setw( 12 ) << x;
      }
      cout << '\n';
   }
}

//===

void print( const string &title, const vec &A )
{
   if ( title != "" ) cout << title << '\n';

   for ( auto x : A )
   {
      if ( abs( x ) < NEARZERO ) x = 0.0;
      cout << setw( 12 ) << x << '\n';
   }
}

//===

matrix matmul( const matrix &A, const matrix &B )
{
   int rowsA = A.size(),   colsA = A[0].size();
   int rowsB = B.size(),   colsB = B[0].size();
   assert( colsA == rowsB );

   matrix C( rowsA, vec( colsB, 0.0 ) );
   for ( int i = 0; i < rowsA; i++ )
   {
      for ( int j = 0; j < colsB; j++ )
      {
         // scalar product of ith row of A and jth column of B
         for ( int k = 0; k < colsA; k++ ) C[i][j] += A[i][k] * B[k][j];
      }
   }
   return C;
}

//===

vec matmul( const matrix &A, const vec &V )
{
   int rowsA = A.size(),   colsA = A[0].size();
   int rowsV = V.size();
   assert( colsA == rowsV );

   vec C( rowsA, 0.0 );
   for ( int i = 0; i < rowsA; i++ )
   {
      for ( int k = 0; k < colsA; k++ ) C[i] += A[i][k] * V[k];
   }
   return C;
}

//===

matrix identity( int n )
{
   matrix I( n, vec( n, 0.0 ) );
   for ( int i = 0; i < n; i++ ) I[i][i] = 1;
   return I;
}

//===

bool GaussElimination( matrix A, vec B, vec &X )
{
   int n = A.size();

   // Row operations
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivot
      if ( doPivot )
      {
         int r = i;
         double maxA = abs( A[i][i] );
         for ( int k = i + 1; k < n; k++ )
         {
            double val = abs( A[k][i] );
            if ( val > maxA )
            {
               r = k;
               maxA = val;
            }
         }
         if ( r != i )
         {
            for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
            swap( B[i], B[r] );
         }
      }

      // Make upper-triangular
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;

      for ( int r = i + 1; r < n; r++ )
      {
         double multiple = A[r][i] / pivot;
         for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
         B[r] -= multiple * B[i];
      }
   }

   // Check last row
   if ( abs( A[n-1][n-1] ) < SMALL ) return false;

   // Back-substitute
   X = B;
   for ( int i = n - 1; i >= 0; i-- )
   {
      for ( int j = i + 1; j < n; j++ )  X[i] -= A[i][j] * X[j];
      X[i] /= A[i][i];
   }

   return true;
}

//===

bool GaussInversion( matrix A, matrix &B )
{
   int n = A.size();
   B = identity( n );

   // Row operations
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivot
      if ( doPivot )
      {
         int r = i;
         double maxA = abs( A[i][i] );
         for ( int k = i + 1; k < n; k++ )
         {
            double val = abs( A[k][i] );
            if ( val > maxA )
            {
               r = k;
               maxA = val;
            }
         }
         if ( r != i )
         {
            for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
            for ( int j = 0; j < n; j++ ) swap( B[i][j], B[r][j] );
         }
      }

      // Make A upper-triangular
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;

      for ( int r = i + 1; r < n; r++ )
      {
         double multiple = A[r][i] / pivot;
         for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
         for ( int j = 0; j < n; j++ ) B[r][j] -= multiple * B[i][j];
      }
   }

   // Check last row
   if ( abs( A[n-1][n-1] ) < SMALL ) return false;

   // Normalise diagonal
   for ( int i = 0; i < n; i++ )
   {
      double multiple = 1.0 / A[i][i];
      for ( int j = i; j < n; j++ ) A[i][j] *= multiple;
      for ( int j = 0; j < n; j++ ) B[i][j] *= multiple;
   }

   // Make A the identity
   for ( int i = n - 1; i > 0; i-- )
   {
      for ( int r = i - 1; r >= 0; r-- )
      {
         double multiple = A[r][i];
         A[r][i] = 0;
         for ( int j = 0; j < n; j++ ) B[r][j] -= multiple * B[i][j];
      }
   }

   return true;
}

== Gaussian elimination to solve AX = B ===
X:
   0.0740413
     2.18256
     3.47124
     -5.9638
     2.50015
     -1.2642
B:
           1
          38
          48
          48
          49
          41
AX:
           1
          38
          48
          48
          49
          41


== Calculate B = Ainverse ==
A:
           1           1           1           1           1           1
          11           9          14           5          -1          -1
           4          10          12           2           2           7
           3           7          12           3           4           1
           8          11           3          -3           4          11
           7          13          12           6           3           1
B:
    0.224991    0.105117   -0.100711    0.014384   0.0582119  -0.0696118
   -0.433594   -0.069845  0.00104907  -0.0605665   0.0223103    0.171558
  -0.0862805   0.0361581   0.0491899   0.0934375  -0.0205618  -0.0891479
    0.655041  -0.0175545   0.0766989   -0.161417   -0.102203   0.0761627
    0.263527  -0.0512647   -0.166453    0.209885   0.0601002  -0.0206085
    0.376314 -0.00261103    0.140226  -0.0957221  -0.0178576   -0.068353
AB:
           1           0           0           0           0           0
           0           1           0           0           0           0
           0           0           1           0           0           0
           0           0           0           1           0           0
           0           0           0           0           1           0
           0           0           0           0           0           1

[ronenp88]

But I can't get accurate results that way, and with integers I get precisely the answer, for small matrices though
anyway thanks I will look into what you have written above...

[lastchance]

@ronenp88,
If I saw the number
24689 / 12344
I wouldn't be able to gauge its magnitude. However, it differs from 2.0 only in the 6th significant figure. Believe me, I've never used lab equipment capable of measuring to that accuracy.

There are some good uses of rational fractions ... but inverting matrices isn't one of them. (And gaussian elimination wouldn't be my method of choice for large matrices either.)

[ronenp88]

I added function that reduce equations (see GitHub) so the numbers will not get too big, but for this matrix shown above, it still gets stuck before finishing and I'm using break points but can't find where it gets stuck, in which loop, any help?

this is the last result when it stuck:

1466    0       0       0       0       -1493   1996
0       -2932   0       0       0       -7359   2904
0       0       733     0       0       -956    3753
0       0       0       2932    0       3267    -21616
0       0       0       0       42895   0       107244
0       0       0       0       0       -42895  54228

[lastchance]

Well, if you are using gaussian elimination (and I assume you are) it is on the last stage of pivoting to reduce the last column. It has managed to get the first 0 above the diagonal element, but not the next one. Looking at the numbers involved, since you have to multiply some of those together you may overflow an int. Alternatively, if you are trying to reduce it by an HCF found by repeated trial then, for large numbers, this is going to take an eternity, especially if they are coprime. Use the standard HCF (aka GCD) algorithm.

At a rough guess it baulks at the row operation
R3 -> 42895 R3 + 3267 R5
where 'Ri' means row i and i is indexed from 0.

But the equation in the last line looks OK:
-42895X5 = 54228, giving X5 = -1.2642...
which is what my code gave earlier. Similarly, X4 = 107244/42895 = 2.50015... as before. So, apart from overflowing, or possibly simply taking an age to factor down, you are probably getting the right answer.

Do what the rest of the world does and use doubles. Inverting matrices is not a good application for rational fractions. Sooner or later you will overflow ints or long longs, and repeated use of an HCF algorithm will take an eternity.

Your code is unnecessarily large, with about 7 levels (edit: 8 levels) of indentation in places and massive amounts of repetition. You are asking rather a lot of forum users when you aren't even posting it here.

[ronenp88]

It was too large to post here, sorry.
and thanks anyway I guess you are right what you said about the algorithm, maybe I should use doubles instead, thanks

[ronenp88]

void reduce_equation(vector<vector<rational>>& matrix, long long int row) {
	long long int biggest_elem = 0;

	for (int i = 0; i < matrix[row].size(); ++i) {
		if (biggest_elem < std::abs(matrix[row][i].num)) {
			biggest_elem = std::abs(matrix[row][i].num);
		}
	}

	for (; biggest_elem >= 2; biggest_elem--) {
		bool state = true;
		for (int i = 0; i < matrix[row].size(); ++i) {
			if (matrix[row][i].num%biggest_elem == 0) {
				continue;
			}
			else {
				state = false;
				break;
			}
		}
		if (state) {
			for (int j = 0; j < matrix[row].size(); ++j) {
				matrix[row][j].num /= biggest_elem;
			}
		}
	}
}

This is GCD right?


BTW I see now the problem, biggest_elem is huge (750055444) which this is why it takes so long, I guess it doesn't stuck, it just takes eternity you are right..

[lastchance]

No idea! Your code is indecipherable.

For two numbers I just use

int hcf( int a, int b )
{
   int q = 1;

   while ( b > 0 )
   {
      q = b;
      b = a % q;
      a = q;
   }
   return q;
}

[ronenp88]

It wasn't the GCD, I changed the code to this and now the calculation finishes:

long long int gcd(long long int a, long long int b) {
	if (a == 0)return b;
	return gcd(b%a, a);
}

void reduce_equation(vector<vector<rational>>& matrix, long long int row) {
	long long int result = matrix[row][0].num;
	for (int i = 1; i < matrix[row].size(); i++)
		result = gcd(matrix[row][i].num, result);

	for (int i = 0; i < matrix[row].size(); ++i) {
		matrix[row][i].num /= result;
	}	
}

X1 = 3176/42895
X2 = 93621/42895
X3 = 148899/42895
X4 = -255817/42895
X5 = 107244/42895
X6 = -54228/42895