OPR-computation-related linear algebra problem

[Nikhil Bajaj]

New Code, based on what James put up (I just added some disp's so that the results would be more clear. disps are outside of tics and tocs. I did not find any bugs though had to change #'s to %'s.

clc
disp('Loading Data...')
tic
d = load('d.dat');
N = load('N.dat');
toc
Ns = sparse(N);
D = diag(Ns);
Ds = sparse(diag(D)); %This was a bug... maybe it still is!

% Reference Solution 
disp('Reference Solution:')
tic
output1 = Ns\d;
toc


% CG Solution
disp('CG Solution:');
tic
output2 = pcg(Ns,d);
toc

% Diagonal PCG Solution
disp('Diagonal PCG Solution:');
tic
output3 = pcg(Ns,d,[],[],Ds);
toc

% Reverse Cutthill-McKee re-ordering
disp('Re-ordering (Reverse Cutthill-McKee:');
tic
p = symrcm(Ns); % permutation array
Nr = Ns(p,p); % re-ordered problem
toc

% Re-ordered Solve
disp('Re-ordered Solution:');
tic
output4 = Nr\d; %answer is stored in a permuted matrix indexed by 'p'
toc

Output:

Loading Data...
Elapsed time is 3.033846 seconds.
Reference Solution:
Elapsed time is 0.014136 seconds.
CG Solution:
pcg stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 4.8e-05.
Elapsed time is 0.007545 seconds.
Diagonal PCG Solution:
pcg converged at iteration 17 to a solution with relative residual 8.9e-07.
Elapsed time is 0.009216 seconds.
Re-ordering (Reverse Cutthill-McKee:
Elapsed time is 0.004523 seconds.
Re-ordered Solution:
Elapsed time is 0.015021 seconds.

I didn't precondition earlier because I was being sloppy/lazy . Thanks for calling me out. And you're right, I should have used pcg. Thanks for the suggestion.

[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

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

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

[Foster]

Originally Posted by RyanCahoon

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

Long, long ago I was the support for a Pascal compiler that ran on the Unisys A-Series. That compiler did Row major, since that was the underlying architecture of a stack machine. A 2D array consisted of an array of pointers, each pointer referencing a row. Once you got the row, it was fast to "walk the row" since the elements were contiguous in memory.

Walking the matrix in col order caused two memory accesses per read, one for the pointer (and the time to do the reference fetch) and the read of the data elements. You could see a -10x reduction in speed doing column order vs row order. -- Thanks for the "back when" reminder.

Ether: Just for grins, I tried RLab on our 16 way cluster. For some reason I'm not getting responses to tic()/toc(); What Windows OS are you running RLab under?

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

[Joe Ross]

I just got a new computer at work with a Xeon E5-1620 and a Quadro 4000 GPU (256 CUDA cores). Using Matlab 2012b:

CPU Invert and Multiply: 1.4292s
CPU Linear Solver: 0.22521s
CPU Sparse Linear Solver: 0.034423s

GPU Invert and Multiply: 0.537926s
GPU Linear Solver: 0.218403s

Loading the N matrix into the GPU was 1.890602s.
Creating the Sparse Matrix for N was 0.032979s.

[RyanCahoon]

Originally Posted by Ether

Ryan, could you please re-run this, without iterating? I want to eliminate the possibility that Matlab is optimizing out the iteration.

I had tried this originally, and the results were consistent with the iterated/averaged result, but I was getting some variation in timing so I wanted to take the average case. Interestingly, the average from the iterated trials was consistently higher than any of the trials running single-shot.

Linear solver 0.19269
Invert and multiply 1.8698

N = dlmread('N.dat');
d = dlmread('d.dat');

tic;
    r1 = N \ d;
t1 = toc;

% also save r1 to a file here so the computation is not optimized out.
dlmwrite('r_solver.dat', r1);

disp(['Linear solver ' num2str(t1)]);


tic;
    r2 = inv(N) * d;
t2 = toc;

% also save r2 to a file here so the computation is not optimized out.
dlmwrite('r_invmult.dat', r2);

disp(['Invert and multiply ' num2str(t2)]);

[flameout]

Originally Posted by Ether

PS - can someone with a working Octave installation please run this? also SciLab and R

Since no-one has done Octave yet, I'll go ahead and do it (along with MATLAB for comparison). I can't do SciLab or R because I don't know how to use those

MATLAB 2012b:

>> N = dlmread('N.dat');
>> d = dlmread('d.dat');
>> tic ; r = N \ d; toc
Elapsed time is 0.797772 seconds.

GNU Octave 3.6.2:

octave:1> N = dlmread('N.dat');
octave:2> d = dlmread('d.dat');
octave:3> tic ; r = N \ d; toc
Elapsed time is 0.624047 seconds.

This is on an Intel i5 (2 core + hyperthreading) with Linux as the host OS (kernel version 3.7.6).

[Ether]

Originally Posted by flameout

Since no-one has done Octave yet, I'll go ahead and do it

Thanks flameout.

I can't do SciLab or R because I don't know how to use those

Just installed SciLab 5.4.1 with Intel Math Kernel Library 10.3 on a 7-year-old desktop PC:

• Intel Pentium D 3.4GHz (x86 Family 15 Model 6 Stepping 4)
  
• 32-bit XP Pro SP3
  
• 500GB Seagate Barracuda 7200

-->stacksize(70000000);
 
-->tic; N=read("N.dat",2509,2509); toc
 ans  = 1.672  
 
-->d=read("d.dat",2509,1);
 
-->tic; x=N\d; toc
 ans  = 1.672  
 
-->tic; Ns=sparse(N); toc
 ans  = 0.141  
 
-->tic(); xs = umfpack(Ns,'\',d); toc
 ans  = 0.14

[flameout]

Originally Posted by Ether

Just installed SciLab 5.4.1 with Intel Math Kernel Library 10.3 on a 7-year-old desktop PC:

Now that I have working SciLab code, I'll go ahead and re-do the tests (with additional instrumentation on the reading). This is on the same computer as before (dual-core, hyperthreaded Intel i5, Linux 3.7.6).

MATLAB R2012b:

>> tic ; N = dlmread('N.dat'); toc
Elapsed time is 3.074810 seconds.
>> tic ; d = dlmread('d.dat'); toc
Elapsed time is 0.006744 seconds.
>> tic ; r = N \ d; toc
Elapsed time is 0.323021 seconds.
>> tic ; dlmwrite('out.dat', r); toc
Elapsed time is 0.124947 seconds.

I noticed that the solve time was very different from my previous run. This may be due to dynamic frequency scaling -- the results in this post (for all software) is with the frequency locked at the highest setting, 2.501 Ghz. It may also be due to a disk read -- I had not loaded MATLAB prior to running the previous test; now it's definitely in the disk cache. The solve is now consistently taking the time reported above, about a third of a second.

GNU Octave 3.6.2:

octave:1> tic ; N = dlmread('N.dat'); toc
Elapsed time is 1.87042 seconds.
octave:2> tic ; d = dlmread('d.dat'); toc
Elapsed time is 0.00241804 seconds.
octave:3> tic ; r = N \ d; toc
Elapsed time is 0.528489 seconds.
octave:4> tic ; dlmwrite('out.dat', r); toc
Elapsed time is 0.00820613 seconds.

Octave seems more consistent. The solve time is higher than for MATLAB, but the I/O times are consistently better.

Scilab 5.3.3:

-->stacksize(70000000);
 
-->tic; N=read("N.dat", 2509, 2509); toc
 ans  =
 
    1.21  
 
-->tic; d=read("d.dat", 2509, 1); toc
 ans  =
 
    0.003  
 
-->tic; x=N\d; toc
 ans  =
 
    1.052  
 
-->tic; Ns=sparse(N); toc
 ans  =
 
    0.081  
 
-->tic(); xs = umfpack(Ns,'\',d); toc
 ans  =
 
    0.081

Scilab failed to read the provided d.dat out of the box (reporting an EOF before it was done reading 2509 rows). I was able to correct this by adding a single newline to the end of d.dat.

FreeMat 4.0:

--> tic ; N = dlmread('N.dat'); toc
ans =
    2.8630
--> tic ; d = dlmread('d.dat'); toc
ans =
    0.0080
--> tic ; r = N \ d; toc
ans =
    3.4270

FreeMat did not have a dlmwrite function, so I haven't reported the write times for it. The time it took to solve the equations was significantly slower than any of the other programs. This did not improve with subsequent runs.