C --> Timing a Matrix Inversion

[font=courier new][size=8pt]Here is a C program which fills an N x N matrix with random doubles, and then,
inverts it, using the Gauss-Jordan Method.

It performs the inversion, and then calculates and writes an N-dimensional check
vector, to an output file. If the check vector contains all '1's, then, the
inversion is correct.

I compiled it using Bloodshed Dev-C++, version 4.9.9.2.

http://www.bloodshed.net/

The compiled '.exe' file, is only 20 KB.

On my computer, which was manufactured in 1945 (maybe a slight exaggeration), it
takes in the range of, 33 to 37 seconds, to fill a 1000 x 1000 matrix with
random doubles, and to invert it.

I attached the source file below, in case anyone wants to try running the
program. I called the attachment file, "main.txt". If you download it, then,
change its name to "main.c". Put "main.c" into some folder. In the Bloodshed
IDE, you have to make a new project (File --> New --> Project). You can call the
project, anything you want. In the dialog box, choose "C", not, "C++"; and,
"Console application". Then, choose the folder that contains "main.c", as the
project folder. The IDE will then create and open a default file, called,
"main.c". Close that file, without saving it. Then, add your version of "main.c"
(from inside the project folder), to the project (Project --> Add to Project).
Then, open, "main.c". Now, to run the program, just choose, "Execute --> Compile
& Run" (F9). When you run it, a console window will open. The program will write
its progress in the window, and will tell you when it finishes. After it
finishes running, you can open the output file (by default, called,
"matinv.txt"), in the IDE. That file gives the run's statistics (matrix size,
elapsed time), and lists the calculated N-dimensional check vector. If the
inversion is correct, then all of the N check vector values, should be very
close to 1.

:twisted:
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 setpivot(double m1[N][N], double m2[N][N], int col);
void switchrows(double m1[N][N], int row1, int row2);
void dividerow(double m1[N][N], int row, double divisor);
void torowaddrow(double m1[N][N], int torow, int addrow, double addrowfactor);
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 row, col;
double divisor, factor;

for(col = first; col <= last; col++)
{
setpivot(original, inversion, col);
divisor = original[col][col];
dividerow(original, col, divisor);
dividerow(inversion, col, divisor);

for(row = col + 1; row <= last; row++)
{
factor = -original[row][col];
torowaddrow(original, row, col, factor);
torowaddrow(inversion, row, col, factor);
}
}
for(col = last; col >= 1; col--)
for(row = col - 1; row >= first; row--)
{
factor = -original[row][col];
torowaddrow(original, row, col, factor);
torowaddrow(inversion, row, 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(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 setpivot(double m1[N][N], double m2[N][N], int col)
{
int row, maxrow;
double max = 0;
double test;

for(row = col; row <= last; row++)
{
test = abs(m1[row][col]);
if( test > max)
{
max = test;
maxrow = row;
}
}
if(maxrow != col)
{
switchrows(m1, col, maxrow);
switchrows(m2, col, maxrow);
}
}

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

void switchrows(double m1[N][N], int row1, int row2)
{
double temp;
int i;
for(i = first; i <= last; i++)
{
temp = m1[row1][i];
m1[row1][i] = m1[row2][i];
m1[row2][i] = temp;
}
}

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

void dividerow(double m1[N][N], int row, double divisor)
{
int i;
for(i = first; i <= last; i++) m1[row][i] /= divisor;
}

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

void torowaddrow(double m1[N][N], int torow, int addrow, double addrowfactor)
{
int i;
for(i = first; i <= last; i++) m1[torow][i] += m1[addrow][i] * addrowfactor;
}

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

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()
{
time_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].\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]

dev-c++ is very nice. If you want to make GUI c++ apps you might want to look at wxDev-c++
It is dev-c++ with wxWidgets packaged in with nice wizards to get you started.
http://wxdsgn.sourceforge.net/

running Windows XP(32-bit) under VMware Fusion on my mac pro it takes about 21 seconds, the same program compiles without complaint on the mac and runs in about 15 seconds but the check vector is always full of NaN's

[font=courier new][size=8pt]I don't know where the NaNs are coming from, jack. When I run it, the check vector is always filled with values, each of which is very very close to 1.

From what you wrote, I am guessing that the check vector was correct when you ran it on your Mac using XP, is that true? I hope so.

But, when you compiled it with a Mac C compiler, you got the bad check vectors? If so, I don't know why that would be.

:oops: :(
Dan

(I see that wxDev-C++ is still being developed, Kent; --> that's good. But, it seems that nothing has been done with Dev-C++, since 2005; --> maybe, that's not so good. On the other hand, I guess C and C++ have not changed much lately.)

Dan, yes to both of your questions, I will have to take a look at code.

Dan, the problem on the mac is that the random numbers were 32-bit instead of 16.

[font=courier new][size=8pt]I put in the code for the output file, that the random numbers vary from 0 to 32767. I should have added that this range was for my machine (I checked RAND_MAX for my machine.), and not necessarily for every machine. According to "C in a Nutshell", RAND_MAX is equal to at least 32767, but, it can be bigger. It varies according to the machine architecture, I guess. But, the matrix that the random numbers are assigned to, is of type double (64 bits). So, whether the random numbers are 16 bits or 32 bits, I think the compiler should automatically convert them to doubles. Unless, for your Mac C compiler, you have to provide an explicit cast in the code. In that case, I guess the line in the function initializeall(), would change from,

m1{i}{j} = rand();

to,

m1{i}{j} = (double) rand();

(I substituted brackets ({}) for braces (&#91;]), above. The output of the forum editor doesn't display braces properly, probably because they are used to format text.)

Concerning the explicit cast, I could be wrong, I'm no expert.

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

Anyway, I'm glad that you at least got the program to work on your Mac, using the virtual XP environment.

:x :P
Dan

[font=courier new][size=8pt]I was able to reduce the time it takes to fill a 1000 x 1000 matrix with random
doubles, and to invert it. On my machine, the reduction is approximately 26%
(previously, the operation averaged 35 seconds on my machine, and now, it
averages 27 seconds.)

:( :P :x
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 = abs(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(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()
{
time_t ticks;
double et;
ticks = clock();
et = (double)ticks / CLOCKS_PER_SEC;
return et;
}

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

void writeoutfile(double t, double v1[N])
{
int i, j;
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]

You guys got some fast PC's. I get 30.14 secs for the first program and 29.56 for the second.

Dan, I compiled your first example with gcc with the command
gcc main.c -o main
execution time on the mac was about 15 seconds, with optimization it was about 7 seconds, under windows the speedup was only a few seconds.
gcc main.c -O3 -o main

under Windows XP(32) your second example runs in about 18 second unoptimized and about 13 seconds optimized
on the Mac it runs in about 12 and 7 seconds repectivly.

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

We are left in the dust.

:oops: :( :unguee:
Dan

[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.)

:oops: :x :twisted:
Dan

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!

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.

[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.

:oops: :x :twisted:
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]

[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

:oops: :x :twisted:
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]

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)