.NET Opensource

I have an application that uses an equation solver.  The existing version is written in VB6 and it seems to work quite well but it is a little slower than I would like.  I found Math.NET via a web search and was very pleased to see that it had a NUnit test suite with it.  I have extend that test suite to test the Solve method of the Matrix class.  Some of the tests are small hard coded matrices and some are dynamically created filled with random numbers.  My intent was to find out if the solution time was any better than my own solver.  Unfortunately my quick and dirty test doesn't behave very well.  What I do is create an NxN matrix of random doubles and an Nx1 vector of ranom doubles and multiply them to get the right hand side then I use the original square matrix and the new right had side in the Solve method to get a vector that should be the same as the original 1xN  vector.

It looks like it works for small matrices but larger ones fail seemingly randomly.  I'm willing to believe that there are good reasons for the failure, including my lack of understanding of the system.  The problem is that the Solve method does not throw any exceptions, it just silently gives the wrong answer.

So at last to my first question: are there any known limits to the size of matrix that the Solve method can handle?

Secondly are they any ways of detecting failures? 

 I need to handle matrices of at least 330x330 elements.  These represent thermal networks so many of the entries are zero which means that my naive tests are not representative of the real problems I want to solve.

(Of course now I wonder how well my own code really performs, have to write a parallel test suite) 

 Just in case they might be of use to you I have attached the two classes that I added to the test suite.

Thank you very much for the report.

Indeed, I was able to reproduce the problem. It usually started failing at TestMatrix_Solve080.

However, luckily I was able to find the problem: In TestMatrix you didn't pass forward the parameter epsilon to the sub-procedure TestMatrix_Solutions, but always used 1e-13 instead. Apparently that precision was not reached any more starting from about 80x80, hence the tests started to fail. I changed TestMatrix to pass on the epsilon parameter, and all tests went green.

Btw, even though the precision doesn't reach 1e-13, a test with a 1000x1000 matrix usually worked for me with precision 1e-9, and required 2145 ms on my machine (debug code, through VS->TestDriven->NUnit):

TestMatrix_NxN(1000, 1e-9);

Back to the questions:

a) At least I'm not aware of any specific. There are of course: at least one is given by the Int32 limits as we use integers to index the arrays. I would have to scan through the code base to find out the exact theoretical limits. I just tried to solve a 10000x10000 matrix but got back an out of memory exception. That would be a 200 MB array if I calculated correctly, should theoretically be possible, but apparently the .Net CLR thinks different.

b) Failures should throw exceptions. If you want to detect overflows (that normally do not throw exceptions in .Net) you may have to use the checked-keyword however. (Update: I think checked is only relevant for integers, as floating point numbers overflow to infinity anyway - but I might be wrong on this)

Btw, you mentioned that many of your values are zeros. Iridium actually once had a sparse-matrix implementation but I temporarily removed due to some integration problems (and lack of time to fix them).

Update: In case you're interested, I usually get the following times here (again with debug code):

100x100: 3 ms
200x200: 20 ms
330x330: 100 ms
1000x1000: 2200 ms

Christoph Rüegg:

However, luckily I was able to find the problem: In TestMatrix you didn't pass forward the parameter epsilon to the sub-procedure
TestMatrix_Solutions, but always used 1e-13 instead.

... 

Update: In case you're interested, I usually get the following times here (again with debug code):

100x100: 3 ms
200x200: 20 ms
330x330: 100 ms
1000x1000: 2200 ms

Oops!  That is so embarrassing.  Thanks for taking the time to find my bug.  And thanks for the timings too. 

Looks like you have a faster machine than me (1.8GHz T32p IBM laptop).  I get this:

10x10               4
60x60                5
100x100           11
300x300         147
400x400         324
600x600       1179
800x800       2666
1000x1000    5230

Not the dramatically better timing (in comparison with my own VB6) I was hoping for but enough to justify further investigation; in particular solving a realistic problem. 

Red facedly yours,

Kevin Whitefoot. 

Yes, it's a 3.4 GHz P4 desktop

I've just tested it also on my notebook (2GHz):

Debug Code:

100x100: 5 ms
200x200: 40 ms
330x330: 160 ms
1000x1000: 4000 ms

Release Code:

100x100: 1 ms
200x200: 20 ms
330x330: 70 ms
1000x1000: 1900-1950 ms

Apparently the code built in release mode is about twice as fast as the debug code.

It might be a good idea to take advantage of Microsoft's Accelerator Project that takes use of the GPU. This should speed up the numeric matrix calculations.

Download: 
http://research.microsoft.com/research/downloads/Details/50ee362a-c4d7-4fe6-9018-1b7f9c1dd5dc/Details.aspx

Wiki:
http://channel9.msdn.com/wiki/default.aspx/Accelerator.HomePage 

Looks like making a stupid mistake can be good sometimes.  I would never have posted here if my little benchmarks had worked and I wouldn't have  got the very useful bits of information  I have.

 I'll try again with release instead of debug and the GPU idea sounds very interesting.  Unfortunately the link on the MS page gives a 550 file not found error!  This link worked though:

http://research.microsoft.com/research/downloads/Details/25e1bea3-142e-4694-bde5-f0d44f9d8709/Details.aspx?CategoryID.

Thanks.
 

Well, there keeps the problem that Accelerator needs DirectX 9 which is only available on Windows but not on Unix based systems.

So, even you can use Accelerator with Mono you can't use it on Linux or a Mac which makes this a insular solution and therefore requires double coding - one time for Windows and one time for Unix.