C --> Timing a Matrix Inversion

[danbaron]

[font=courier new][size=8pt]jack has the fast machines, Kent.

We are left in the dust.


Dan

[danbaron]

[font=courier new][size=8pt]I didn't even know that I had "gcc" on my (Windows) machine. I don't remember when or how I got it.

Anyway, I was unsure if maybe I did have it, so, from the the command line I typed, "gcc", and it responded, "gcc: no input files". So, then I knew I had it.

So, I compiled both versions, using jack's command,

"gcc main.c -O3 -o main".

I guess that "O3", must mean, super optimization.

I ran both versions three times. Here are the results (in seconds):

slow version --> 16.076, 15.334, 15.084
fast version --> 11.877, 11.744, 11.872

gcc is good, and, free.

(Now I use your tip, Kent, to open a command window in the folder where the source files reside (shift - right click). Previously, I would open a command window, and then type, "cd c:\users\root\desktop\c", for instance. And, many times, I would have a spelling error, and then have to re-type most of the line - very primitive, very frustrating.)


Dan

[kryton9]

I will need to play with those commands too. I never used them. I always use the IDE's to compile c++, and set the compilation flags that way. I usually set it for optimal size instead of speed.
Just recompiled using that great parameter that Jack and Dan mentioned, wow.

First one 8.961, second test 7.630 filesize about 32K. My compiled for optimal size, but very slow speed filesize about 9K.

My old dog computer gets a new bark after all, that makes me very happy. Thanks guys!

[MikeStefanik]

For fun, I compiled this using Visual Studio 2010 Profesisonal (release build, which is optimized), elapsed time was 2.652 seconds.

Edit: By the way, you have an unreferenced variable (j) in main, your code is assuming that time() returns a 32-bit value (not always the case) and in C you should use fabs() instead of abs() since it's working with a double (abs() is for ints). Not a big deal if you're using C++, since abs() is typically overloaded, but just so you know.

[danbaron]

[font=courier new][size=8pt]I am the same way, Kent, I have gotten used to using IDEs. It's good that jack
reminded us about the power we have been neglecting due to our laziness.

------------------------------------------

I think I found the extraneous "j", Mike. It was in, writeoutfile().

------------------------------------------

The code where time() appears, is,

srand(time(&t));

In, "C in a Nutshell", it says that srand(), takes an argument of type,
"unsigned". But, it doesn't say which one, I think there are four, "unsigned
short", "unsigned int", "unsigned long", and, "unsigned long long" (<-- C99).
time() returns the type, time_t, "generally as long or unsigned long". So, I'll
change the code to this,

srand((unsigned long)time(&t));

Presumably, "unsigned long", corresponds to 4 bytes, yes?.
Anyway, I don't know what else to do.

------------------------------------------

I looked it up, and you are absolutely correct, about using fabs() instead of
abs().

------------------------------------------

The 2.652 seconds, using Visual Studio 2010, is amazingly fast.

However, believe it or not, it is not the current champion.

If you look here (next link), at this little bit of Python code, you will see that the
timing is for a 6000 x 6000 matrix.

http://community.thinbasic.com/index.php?topic=3485.0

I changed the code to for a 1000 x 1000 matrix, and ran it three times from the
command line. The times were (seconds),

1.730
1.725
1.714

------------------------------------------

I will post the C code below again, with the changes due to Mike.


Dan

[code=c]//------------------------------------------------------------------------------------------------

// A 'C' program to time matrix inversions.
// Uses the Gauss-Jordan Method to invert an N x N matrix, filled with random doubles.
// The results are written to a file; by default, "matinv.txt".

//------------------------------------------------------------------------------------------------

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//------------------------------------------------------------------------------------------------

// Set size of N x N matrix to be inverted.
#define N 1000

//------------------------------------------------------------------------------------------------

const int first = 0;
const int last = N - 1;

// Name of outfile.
const char outfilename[80] = "matinv.txt";

//------------------------------------------------------------------------------------------------

void invert(double original[N][N], double inversion[N][N]);
void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N]);
void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N]);
void matvecmult(double m1[N][N], double v1[N], double v2[N]);
double getelapsedtime();
void writeoutfile(double t, double v1[N]);

//------------------------------------------------------------------------------------------------

// Matrix to be inverted (filled with random doubles).
double a0[N][N];

// Copy of matrix to be inverted.
double am[N][N];

// Calculated inversion of a0.
double ai[N][N];

// Calculated identity matrix.
double id[N][N];

// va is multiplied by id, to give the check vector, vb.
double va[N];

// check vector.
double vb[N];

//------------------------------------------------------------------------------------------------

int main(int argc, char *argv&#91;])
{
double t0, t1, tt;

printf("Inverting matrix..\n");
printf("\n");

t0 = getelapsedtime();

// a0 is filled with random doubles.
// am is set equal to a0.
// ai is set to the N x N identity matrix.
// va is set to an N-dimensional vector, of all '1's.
initializeall(a0, am, ai, va);

// am becomes the N x N identity matrix.
// ai becomes the (N x N) inversion of a0.
invert(am, ai);

t1 = getelapsedtime();
tt = t1 - t0;

printf("Finished inverting a %d by %d matrix.\n", N, N);
printf("elapsed time = %6.3f seconds\n", tt);
printf("\n");
printf("Calculating check vector..\n");
printf("\n");

// Multiply the inverted matrix, ai, by the original matrix, a0, to give the calculated identity matrix, id.
matmatmult(a0, ai, id);

// Multiply the N-dimensional vector of all '1's, va, by the calculated N x N identity matrix, id, to give the check vector, vb.
// If the inversion is correct, then, vb, should also be an N-dimensional vector, of all '1's.
// vb, is written to the output file.
matvecmult(id, va, vb);

writeoutfile(tt, vb);

printf("Finished calculating check vector.\n");
printf("Output written to file, \"%s\".\n", outfilename);
printf("\n");
printf("Press 'Enter' to quit.\n");
getchar();
return 0;
}

//------------------------------------------------------------------------------------------------

void invert(double original[N][N], double inversion[N][N])
// Uses the Gauss-Jordan Method.
{
int pivot, maxrow, row, col;
double max, test, temp, divisor, factor;
//---------------------------------------------
for(pivot = first; pivot <= last; pivot++)
{
//---------------------------------------------
for(row = pivot; row <= last; row++)
{
max = 0;
test = fabs(original[row][pivot]);
if( test > max)
{
max = test;
maxrow = row;
}
}
//---------------------------------------------
if(maxrow != pivot)
{
for(col = first; col <= last; col++)
{
temp = original[pivot][col];
original[pivot][col] = original[maxrow][col];
original[maxrow][col] = temp;
temp = inversion[pivot][col];
inversion[pivot][col] = inversion[maxrow][col];
inversion[maxrow][col] = temp;
}
}
//---------------------------------------------
divisor = original[pivot][pivot];
for(col = first; col <= last; col++)
{
original[pivot][col] /= divisor;
inversion[pivot][col] /= divisor;
}
//---------------------------------------------
for(row = first; row <= last; row++)
{
if(row != pivot)
{
factor = -original[row][pivot];
for(col = first; col <= last; col++)
{
original[row][col] += original[pivot][col] * factor;
inversion[row][col] += inversion[pivot][col] * factor;
}
}
}
}
}

//------------------------------------------------------------------------------------------------

void initializeall(double m1[N][N], double m2[N][N], double m3[N][N], double v1[N])
{
int i, j;
time_t t;

srand((unsigned long)time(&t));

for(i = first; i <= last; i++)
for(j = first; j <= last; j++)
{
m1[i][j] = rand();
m2[i][j] = m1[i][j];
}

for(i = first; i <= last; i++)
{
v1[i] = 1;
for(j = first; j <= last; j++)
{
m3[i][j] = 0;
if(i == j) m3[i][j] = 1;
}
}
}

//------------------------------------------------------------------------------------------------

void matmatmult(double m1[N][N], double m2[N][N], double m3[N][N])
{
int i, j, k;
for(i = 0; i <= last; i++)
{
for(j = 0; j <= last; j++)
{
m3[i][j] = 0;
for(k = 0; k <= last; k++) m3[i][j] += m1[i][k] * m2[k][j];
}
}
}

//------------------------------------------------------------------------------------------------

void matvecmult(double m1[N][N], double v1[N], double v2[N])
{
int i, j;
for(i = 0; i <= last; i++)
{
v2[i] = 0;
for(j = 0; j <= last; j++) v2[i] += m1[i][j] * v1[j];
}
}

//------------------------------------------------------------------------------------------------

double getelapsedtime()
{
clock_t ticks;
double et;
ticks = clock();
et = (double)ticks / CLOCKS_PER_SEC;
return et;
}

//------------------------------------------------------------------------------------------------

void writeoutfile(double t, double v1[N])
{
int i;
FILE *outfile;

outfile = fopen(outfilename, "w");
fprintf(outfile, "Inversion of a %d by %d matrix,\n", N, N);
fprintf(outfile, "filled with uniformly distributed random doubles,\n");
fprintf(outfile, "in the range of, [0, 32767] (on my machine).\n");
fprintf(outfile, "\n");
fprintf(outfile, "elapsed time = %6.3f seconds\n", t);
fprintf(outfile, "\n");
fprintf(outfile, "If the inversion was calculated correctly, then every\n");
fprintf(outfile, "value shown below should be very close to 1.\n");
fprintf(outfile, "\n");

for(i = first; i <= last; i++) fprintf(outfile, "%5i %19.16f\n", i + 1, v1[i]);

fclose(outfile);
}

//------------------------------------------------------------------------------------------------
[/code]

[danbaron]

[font=courier new][size=8pt]Here is a new C version.

It uses the same algorithm that I used for the PowerBASIC version. I don't think
I can improve the algorithm.

By the way, there were errors in the previous C versions. They resulted in the
subroutine not performing pivoting. You can invert a matrix without pivoting,
but, if the matrix is close to being singular, or, if there is a zero on the main
diagonal, then the inversion will be incorrect. Since, the matrices here were filled
with random values, there was little danger of that, and the results were OK. I'll
show the part of the subroutine with the errors now.

[code=c]void invert(double original[N][N], double inversion[N][N])
// Uses the Gauss-Jordan Method.
{
int pivot, maxrow, row, col;
double max, test, temp, divisor, factor;
//---------------------------------------------
for(pivot = first; pivot <= last; pivot++)
{
//---------------------------------------------
for(row = pivot; row <= last; row++)
{
max = 0;
test = fabs(original[row][pivot]);
if( test > max)
{
max = test;
maxrow = row;
}
}
[/code]

[font=courier new][size=8pt]Instead, the code part should have been this.

[code=c]void invert(double original[N][N], double inversion[N][N])
// Uses the Gauss-Jordan Method.
{
int pivot, maxrow, row, col;
double max, test, temp, divisor, factor;
//---------------------------------------------
for(pivot = first; pivot <= last; pivot++)
{
max = 0;
maxrow = pivot;
//---------------------------------------------
for(row = pivot; row <= last; row++)
{
test = fabs(original[row][pivot]);
if( test > max)
{
max = test;
maxrow = row;
}
}
[/code]

[font=courier new][size=8pt]Anyway, the previous versions don't matter. The new version, shown below, is
better. I attached the file below, and called it, "matinv.txt". If you download it,
then change the name to "matinv.c", before you compile it.

If you have, "gcc", then, you can compile and run the program, like this.

Open a command window in the folder where the file resides ("shift - right
click" on folder).

Compile the program with the command,
"gcc matinv.c -O3 -o matinv".

Run it with the command,
"matinv".

I ran it three times on my machine. Here is the comparison with the previous
(fast) version in seconds,

previous new
11.877 7.905
11.744 7.876
11.872 7.774


Dan

[code=c]//************************************************************************************************

// FILE = "matinv.c"

// A 'C' program to time matrix inversions.
// Uses the Gauss-Jordan Method to invert an N x N matrix, filled with random doubles.
// The results are written to a file; by default, "matinv.txt".

//************************************************************************************************

// If you have, "gcc", then, you can compile and run the program, like this.
// Open a command window in the folder where this file resides ("shift - right
// click" on folder).

// Compile the program with the command,
// "gcc matinv.c -O3 -o matinv".

// Run it with the command,
// "matinv".

//************************************************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//************************************************************************************************

// Set size of N x N matrix to be inverted.
#define N 1000

//************************************************************************************************

// Name of outfile.
const char outfilename[80] = "matinv.txt";

//************************************************************************************************

void invert(int n, double a[n][n]);
void initializeall(int n, double m1[n][n], double m2[n][n], double v1[n]);
void matmatmult(int n, double m1[n][n], double m2[n][n], double m3[n][n]);
void matvecmult(int n, double m1[n][n], double v1[n], double v2[n]);
double getelapsedtime();
void writeoutfile(int n, double t, double v1[n]);

//************************************************************************************************

// My experience is that, you have to declare these matrices and vectors globally. Otherwise, for
// big values of N (above), Windows will say that the program has stopped working. I guess,
// the function stack is too small. I don't know how to make it bigger.

// Matrix to be inverted (filled with random doubles).
double a0[N][N];

// Calculated inversion of a0.
double ai[N][N];

// Calculated identity matrix.
double id[N][N];

// va is multiplied by id, to give the check vector, vb.
double va[N];

// check vector.
double vb[N];

//************************************************************************************************
//************************************************************************************************

int main(int argc, char *argv&#91;])
{
const int n = N;
double t0, t1, t2, tt;

//---------------------------------------------

printf("\n");
printf("\n");
printf("Inverting a %d by %d matrix..\n", n, n);
printf("\n");

// a0 is filled with random doubles.
// ai is set equal to a0.
// va is set to an N-dimensional vector, of all '1's.
initializeall(n, a0, ai, va);

t0 = getelapsedtime();
// ai becomes the (N x N) inversion of a0.
invert(n, ai);
t1 = getelapsedtime();
tt = t1 - t0;

printf("Finished.\n", N, N);
printf("\n");
printf("elapsed time = %6.3f seconds\n", tt);
printf("\n");
printf("Calculating check vector..\n");
printf("\n");

// Multiply the inverted matrix, ai, by the original matrix, a0, to give the calculated identity matrix, id.
matmatmult(n, a0, ai, id);
// Multiply the N-dimensional vector of all '1's, va, by the calculated N x N identity matrix, id, to give the check vector, vb.
// If the inversion is correct, then, vb, should also be an N-dimensional vector, of all '1's.
// vb, is written to the output file.
matvecmult(n, id, va, vb);

writeoutfile(n, tt, vb);

printf("Finished.\n");
printf("\n");
printf("Output written to file, \"%s\".\n", outfilename);
printf("\n");
printf("Press \"Enter\" to quit.\n");
getchar();
return 0;
}

//************************************************************************************************
//************************************************************************************************

void invert(int n, double a[n][n])
// Uses the Gauss-Jordan Method.
// Inverts matrix a&#91;]&#91;], "in-place".

{
int pivot, row, col, maxrow;
double test, colmax, divisor, factor, temp;
int switchcount, rowswitch, rowswitches[N][2];
int first, last;

//---------------------------------------------

first = 0;
last = n - 1;
switchcount = -1;

for(pivot = first; pivot <= last; pivot++)
{
colmax = 0;
maxrow = pivot;

//---------------------------------------------

for(row = pivot; row <= last; row++)
{
test = fabs(a[row][pivot]);

if( test > colmax)
{
colmax = test;
maxrow = row;
}

}

//---------------------------------------------

if(maxrow != pivot)
{

for(col = first; col <= last; col++)
{
temp = a[pivot][col];
a[pivot][col] = a[maxrow][col];
a[maxrow][col] = temp;
}

++switchcount;
rowswitches[switchcount][0] = pivot;
rowswitches[switchcount][1] = maxrow;

}

//---------------------------------------------

divisor = a[pivot][pivot];
a[pivot][pivot] = 1;

for(col = first; col <= last; col++)
{
a[pivot][col] /= divisor;
}

//---------------------------------------------

for(row = first; row <= last; row++)
{
if(row != pivot)
{
factor = -a[row][pivot];
a[row][pivot] = 0;
for(col = first; col <= last; col++)
{
a[row][col] += a[pivot][col] * factor;
}
}
}

}

//---------------------------------------------

for(rowswitch = switchcount; rowswitch >= 0; rowswitch--)
for(row = first; row <= last; row++)
{
temp = a[row][rowswitches[rowswitch][0]];
a[row][rowswitches[rowswitch][0]] = a[row][rowswitches[rowswitch][1]];
a[row][rowswitches[rowswitch][1]] = temp;
}

}

//************************************************************************************************
//************************************************************************************************

void initializeall(int n, double m1[n][n], double m2[n][n], double v1[n])
{
int i, j;
time_t t;
double d;
int first, last;

first = 0;
last = n - 1;

srand((unsigned long)time(&t));

for(i = first; i <= last; i++)
{
v1[i] = 1;
for(j = first; j <= last; j++)
{
d = rand();
m1[i][j] = d;
m2[i][j] = d;
}
}
}

//************************************************************************************************
//************************************************************************************************

void matmatmult(int n, double m1[n][n], double m2[n][n], double m3[n][n])
{
int i, j, k;
int first, last;

first = 0;
last = n - 1;

for(i = 0; i <= last; i++)
{
for(j = 0; j <= last; j++)
{
m3[i][j] = 0;
for(k = 0; k <= last; k++) m3[i][j] += m1[i][k] * m2[k][j];
}
}
}

//************************************************************************************************
//************************************************************************************************

void matvecmult(int n, double m1[n][n], double v1[n], double v2[n])
{
int i, j;
int first, last;

first = 0;
last = n - 1;

for(i = 0; i <= last; i++)
{
v2[i] = 0;
for(j = 0; j <= last; j++) v2[i] += m1[i][j] * v1[j];
}
}

//************************************************************************************************
//************************************************************************************************

double getelapsedtime()
{
clock_t ticks;
double et;
ticks = clock();
et = (double)ticks / CLOCKS_PER_SEC;
return et;
}

//************************************************************************************************
//************************************************************************************************

void writeoutfile(int n, double t, double v1[n])
{
int i;
int first, last;
FILE *outfile;

first = 0;
last = n - 1;

outfile = fopen(outfilename, "w");
fprintf(outfile, "Inversion of a %d by %d matrix,\n", N, N);
fprintf(outfile, "filled with uniformly distributed random doubles,\n");
fprintf(outfile, "in the range of, [0, 32767] (on my machine).\n");
fprintf(outfile, "\n");
fprintf(outfile, "elapsed time = %6.3f seconds\n", t);
fprintf(outfile, "\n");
fprintf(outfile, "If the inversion was calculated correctly, then every\n");
fprintf(outfile, "value shown below should be very close to 1.\n");
fprintf(outfile, "\n");

for(i = first; i <= last; i++) fprintf(outfile, "%5i %21.19f\n", i + 1, v1[i]);

fclose(outfile);
}

//************************************************************************************************
//************************************************************************************************
[/code]

[zak]

Hi dan
i know nothing about this math subject, but i have tried it, i have MinGW installed, and as you said :"gcc matinv.c -O3 -o matinv", the results is:
elapsed time = 6.204 seconds
but running the main.txt(main.c) in the bottom of the previous page:
elapsed time = 11.031 seconds

i have pentium D cpu 3.0 GHz dual core but of the oldest versions( may be 5 years old or more).
if for some reasons the exe's different, here is my 2 exe's if someone want to try.
http://rapidshare.com/files/411145513/matinv.rar
(virus free, checked with avast)