OPR-computation-related linear algebra problem

[MikeE]

Sorry I'm a bit late to the party.

I'm running Octave 3.2.4, so an older version than flameout
Hardware is Dell E6420 laptop (CPU i7-2640M @ 2.8GHz) Win 7 64bit

octave:2> N = load('N.dat');
octave:3> d = load('d.dat');
octave:4> tic; Ns = sparse(N); toc
Elapsed time is 0.136 seconds.
octave:5> tic; r = Ns \ d; toc
Elapsed time is 0.096 seconds.
octave:6> tic; r = N \ d; toc
Elapsed time is 0.921 seconds.

So solving the sparse matrix is around 230ms and direct solution is around 900ms.

[Ether]

Originally Posted by RyanCahoon

EDIT: Changing the order of the summations got me down to 2.68 and changing to in-place computation like your code got me to 2.58. Beyond that, any improvements would seem to be in the way the Pascal compiler is generating code.

Thanks for doing this Ryan.

Interestingly, the Pascal and C++ compilers I used are essentially identical. Only the front ends are different (for the different languages).

Is it possible that the difference in timing is due to the differences in the memory access due to the data structures we used?

[Ether]

Originally Posted by Greg McKaskle

There is a toolkit called "Multicore Analysis and Sparse Matrix Toolkit", and they ran the numbers using that tool as well.

Greg,

What machine & OS was used?

[Ether]

Originally Posted by BornaE

Loading it into my GTX 580 GPU right now to get some values. Will do that with and without the time taken to load it into the GPU memory and back.

Did you ever get a chance to do this?

[RyanCahoon]

Originally Posted by Ether

Ryan,

Could you edit your post to indicate what O/S you are using?

Windows 7 x64 Home Premium

Originally Posted by Ether

Is it possible that the difference in timing is due to the differences in the memory access due to the data structures we used?

We're both using 8-byte floating point numbers. Assuming Pascal stores a multidimensional array in a contiguous memory block, we're basically using the same data structure - the only difference being that I'm only storing the elements below the diagonal and I believe you store all the elements of the matrix (at least, this is the way that I wrote the read/write code). Maybe this would affect caching.

For comparison, I compiled your code using Free Pascal and the Cholesky decomposition ran in 11.9 seconds on my computer.

[Greg McKaskle]

I haven't used Pascal in a long time, but seem to remember it storing 2D arrays with different elements adjacent. It was column-major and C was row-major. The notation isn't important, but accessing adjacent elements in the cache will be far faster than jumping by 20Kb to pickup up the next element.

Greg McKaskle

[BornaE]

I don't have direct access to my desktop at the moment, I was doing that remotely with Logmein however for some reason I lost the connection and have not got it back yet.

I tried it on a GTX555M with only 24 cuda cores. It was 50% slower than my laptop processor(Core i7 2670QM quad core running at 2.2GHz)
I will post here as soon as I get to my desktop.

I was able to get 0.015 seconds using sparse matrices, however GPU processing does not support sparse matrices directly. I doubt that I can get any faster results than that.

Originally Posted by Ether

Did you ever get a chance to do this?

[Ether]

Originally Posted by BornaE

laptop Core i7 2670QM quad core 2.2GHz...

0.015 seconds using sparse matrices

Matlab, Octave, SciLab, Python, or something else?

Linux, Windows XP/7, 32 or 64 ?

[BornaE]

Matlab 2012b

here are the results

Normal Matrices(CPU and GPU(555M))
Using inv(N)*d:
CPU 1.874s
GPU 2.146s
using N\d:
CPU 0.175s
GPU 0.507s

Sparse Matrices(Only CPU)
Using inv(N)*d:
CPU 0.967s
using N\d:
CPU 0.015s

Cannot get sparse matrices into the GPU easily.

The times are only for the solve operation.

Originally Posted by Ether

Matlab, Octave, SciLab, Python, or something else?

Linux, Windows XP/7, 32 or 64 ?

[Ether]

Originally Posted by Greg McKaskle

I haven't used Pascal in a long time, but seem to remember it storing 2D arrays with different elements adjacent. It was column-major and C was row-major. The notation isn't important, but accessing adjacent elements in the cache will be far faster than jumping by 20Kb to pickup up the next element.

I was thinking the extra computation required for this might be the culprit:

#define ELEMENT(M, i,j) (M[(i)*((i)+1)/2+(j)])

[Ether]

Just installed Octave3.6.4_gcc4.6.2_20130408 on this computer.

Results:

GNU Octave, version 3.6.4

Octave was configured for "i686-pc-mingw32".

octave:1> tic; d = load('k:/MyOctave/d.dat'); toc
Elapsed time is 0.0625 seconds

octave:2> tic; N = load('k:/MyOctave/N.dat'); toc
Elapsed time is 90.0469 seconds

octave:3> tic; r = N \ d; toc
Elapsed time is 1.03125 seconds

octave:4> tic; r = N \ d; toc
Elapsed time is 0.953125 seconds

octave:5> tic; Ns = sparse(N); toc
Elapsed time is 0.0625 seconds
octave:6> tic; rs = Ns \ d; toc
Elapsed time is 0.1875 seconds

octave:7> tic; Ns = sparse(N); toc
Elapsed time is 0.046875 seconds
octave:8> tic; rs = Ns \ d; toc
Elapsed time is 0.171875 seconds

Anybody know why Octave takes so long to load N ?

Load N times, all on the same computer:

Delphi....0.6 seconds

C.........1.3 seconds

SciLab....1.7 seconds

RLab......3.7 seconds

Python....6.5 seconds

Octave...90.0 seconds

[Greg McKaskle]

LV times were on this computer

Win 7 Professional 64-bit
Xeon CPU E5-1620 3.6GHz
16G RAM

Greg McKaskle

[Ether]

RLaB is not a contender for fastest speed, but it's definitely the tiniest linear algebra app out there. It weighs in at under 1.5MB.

Makes a wonderful pop-up for that quick calculation, or for high-school students just learning linear algebra.

Very easy to use... and free.

Welcome to RLaB. New users type `help INTRO'
RLaB version 2.1.05 Copyright (C) 1992-97 Ian Searle

This is free software, and you are welcome to redistribute it under
certain conditions; type `help conditions' for details

> tic(); N=readm("N.dat"); toc()
3.66

> d=readm("d.dat");

> tic(); x=N\d; toc()
28.5

> tic(); x=solve(N,d,"S"); toc()
14.5
>

The second solution method (x=solve(N,d,"S")) tells RLaB that the matrix is symmetric, so it uses LAPACK subroutine DSYTRF (Bunch-Kaufman diagonal pivoting) to solve, which cuts the time in half.

[RyanCahoon]

Originally Posted by Greg McKaskle

I haven't used Pascal in a long time, but seem to remember it storing 2D arrays with different elements adjacent. It was column-major and C was row-major.

I tried reversing the indices and it took 59.1 (vs 11.9) seconds. A couple other websites say Pascal is row-major.

Originally Posted by Ether

I was thinking the extra computation required for this might be the culprit:

#define ELEMENT(M, i,j) (M[(i)*((i)+1)/2+(j)])

I was worried about the same thing. In the most recent version of the code I posted, that macro is only used once (not in a loop) to calculate a pointer to the end of the matrix.

[Ether]

Is anybody else running Octave on a machine with 32-bit-XP Pro?

Are you having the same 30 second delay for Octave to load, and 90 seconds to load the 12MB N.dat file?