Java --> Timing a Matrix Inversion

[danbaron]

[font=courier new][size=8pt]Here is a Java 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 values all very close
to 1, then, the inversion is correct.

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

If anyone wants to try the program, then:

1) Make sure Java (>= 1.5) is installed on your computer, and that it's "bin" directory is in your search path.

1) Download the file, and place it in some folder.

2) Change the name of the file to, "matinv.java".

3) Open a command window in the folder.

4) Perform the command, "javac matinv.java".

5) Perform the command, "java matinv". The program will run in the window. When it finishes, it will write the output file to, "matinv.txt" (inside the folder).

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

I ran it three times for a 1000 x 1000 matrix.

The results were (seconds),

24.096
23.977
24.056

:x :P
Dan


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

import java.lang.System.*;
import java.util.Random;
import java.io.*;

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

public class matinv
{
final static int N = 1000;
final static int first = 0;
final static int last = N - 1;

final static String outfilename = "matinv.txt";

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

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

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

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

// Calculated identity matrix.
static double id[][] = new double[N][N];

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

// check vector.
static double vb[] = new double[N];

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

public static void main(String[] args)
{
BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
long t0, t1;
double tt;
String s;

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

t0 = gettime();

// 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 = gettime();
tt = (t1 - t0) / (double)1000;

System.out.printf("Finished inverting a %d by %d matrix.\n", N, N);
System.out.printf("elapsed time = %6.3f seconds\n", tt);
System.out.printf("\n");
System.out.printf("Calculating check vector..\n");
System.out.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);

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

try
{
s = br.readLine();
}
catch(IOException e)
{
System.out.println("Input error.");
}
}

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

static void invert(double original[][], double inversion[][])
// Uses the Gauss-Jordan Method.
{
int pivot, maxrow, row, col;
double max, test, temp, divisor, factor;

maxrow = -1;

//---------------------------------------------
for(pivot = first; pivot <= last; pivot++)
{
//---------------------------------------------
for(row = pivot; row <= last; row++)
{
max = 0;
test = Math.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;
}
}
}
}
}

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

static void initializeall(double m1&#91;]&#91;], double m2&#91;]&#91;], double m3&#91;]&#91;], double v1&#91;])
{
Random r = new Random();

for(int i = first; i <= last; i++)
for(int j = first; j <= last; j++)
{
m1[i][j] = r.nextDouble();
m2[i][j] = m1[i][j];
}

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

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

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

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

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

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

static long gettime()
{
return System.currentTimeMillis();
}

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

static void writeoutfile(double t, double v1&#91;]) //throws IOException
{
PrintWriter outfile;

try
{
outfile = new PrintWriter(outfilename);

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

for(int i = first; i <= last; i++) outfile.printf("%05d %19.16f\n", i + 1, v1[i]);

outfile.close();
}
catch(IOException e)
{
System.out.println("File Error");
}
}

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

}

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

[danbaron]

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

The previous version had the same errors as the previous versions of the corresponding C program
(http://community.thinbasic.com/index.php?topic=3518.msg26450;topicseen#msg26450), the result being that the inversion
was performed without pivoting. This version fixes those errors, and uses the new and faster algorithm that the
PowerBASIC (http://community.thinbasic.com/index.php?topic=3563.msg26398;topicseen#msg26398) and C programs now use.

It performs and times the inversion, and then calculates and writes an N-dimensional check vector along with the elapsed
time, to an output file. If the check vector contains values all very close to 1, then, the inversion is correct.

I attached the file below, and called it "matinv.txt" (the forum will not permit me to attach a file with the extension,
"java").

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

If anyone wants to try the program, then:

1) Make sure Java (>= 1.5) is installed on your computer, and that it's "bin" directory is in your search path.

1) Download the file, and place it in some folder.

2) Change the name of the file to, "matinv.java".

3) Open a command window in the folder ("shift - right click" on the folder, and choose, "Open command window here",
from the menu).

4) Perform the command, "javac matinv.java".

5) Perform the command, "java matinv". The program will run in the window, indicating its progress. When it finishes, it
will write the output file to, "matinv.txt" (inside the folder).

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

I ran it three times for a 1000 x 1000 matrix.

Here are the results (elapsed time for inversions) for my machine, compared to those of the old version, in seconds.

old new
24.096 12.328
23.977 12.566
24.056 12.319

:oops: :x :twisted:
Dan


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

import java.lang.System.*;
import java.util.Random;
import java.io.*;

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

public class matinv
{

// Set the size of the N x N matrix to be inverted.
final static int N = 1000;

final static String outfilename = "matinv.txt";

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

public static void main(String&#91;] args)
{
BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
long t0, t1;
double tt;
String s;

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

// Matrix to be inverted (filled with random doubles).
double a0&#91;]&#91;] = new double[N][N];

// Calculated inversion of a0.
double ai&#91;]&#91;] = new double[N][N];

// Calculated identity matrix.
double id&#91;]&#91;] = new double[N][N];

// va is multiplied by id, to give the check vector, vb.
double va&#91;] = new double[N];

// check vector.
double vb&#91;] = new double[N];

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

System.out.printf("\n");
System.out.printf("Inverting a %d by %d matrix..\n", N, N);
System.out.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(a0, ai, va);

t0 = gettime();
// ai becomes the (N x N) inversion of a0.
invert(ai);
t1 = gettime();
tt = (t1 - t0) / (double)1000;

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

// Multiply the inverted matrix, ai, by the a 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);

System.out.printf("Finished.\n");
System.out.printf("\n");
System.out.printf("Output written to file, \"%s\".\n", outfilename);
System.out.printf("\n");
System.out.printf("Press 'Enter' to quit.\n");

try
{
s = br.readLine();
}
catch(IOException e)
{
System.out.println("Input error.");
}
}

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

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

{
int pivot, row, col, maxrow;
double test, colmax, divisor, factor, temp;
int switchcount, rowswitch, rowswitches&#91;]&#91;];
int first, last;

rowswitches = new int[N][2];

first = 0;
last = a.length - 1;
switchcount = -1;

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

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

for(row = pivot; row <= last; row++)
{
test = Math.abs(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;
}
}

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

static void initializeall(double m1&#91;]&#91;], double m2&#91;]&#91;], double v1&#91;])
{
int first, last;
double d;
Random r = new Random();

first = 0;
last = m1.length - 1;

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

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

static void matmatmult(double m1&#91;]&#91;], double m2&#91;]&#91;], double m3&#91;]&#91;])
{
int first, last;

first = 0;
last = m1.length - 1;


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

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

static void matvecmult(double m1&#91;]&#91;], double v1&#91;], double v2&#91;])
{
int first, last;

first = 0;
last = m1.length - 1;

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

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

static long gettime()
{
return System.currentTimeMillis();
}

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

static void writeoutfile(double t, double v1&#91;]) //throws IOException
{
int first, last;
PrintWriter outfile;

first = 0;
last = v1.length - 1;

try
{
outfile = new PrintWriter(outfilename);

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

for(int i = first; i <= last; i++) outfile.printf("%05d %19.16f\n", i + 1, v1[i]);

outfile.close();
}
catch(IOException e)
{
System.out.println("File Error");
}
}

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

}

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

[jack]

it runs in about 3.5 seconds on my Mac, the new C version runs in about 3.2 seconds :o