
Hei.
I've been facing a problem while trying to solve a Ax=B matricial system on my VB project. I used to program in Fortran and solve that problem through DGESV LaPack routine, but since I'm new in VB programming I cant seem to get a clue about how to do it.
Thanks in advance, any help would be much appreciated.


Coordinator
Nov 27, 2012 at 11:06 PM

Here's a VB sample:
' Let's solve the matrix equation Ax = b
' Define a 3 X 3 matrix
Dim A As New SquareMatrix(3)
A(0, 0) = 1
A(0, 1) = 1
A(0, 2) = 1
A(1, 0) = 2
A(1, 1) = 1
A(1, 2) = 0
A(2, 0) = 0
A(2, 1) = 1
A(2, 2) = 1
' Define a righthandside vector
Dim b As New ColumnVector(0, 0, 1)
' LU decomposition is the fastest solution method
Dim LU = A.LUDecomposition()
Dim x = LU.Solve(b)
' Show the result
For i As Integer = 0 To x.Dimension  1
Console.WriteLine(x(i))
Next
Note that, in contrast to FORTRAN, array indexes in VB (and other .NET langauges) are zerobased. Note also that there are other approaches besides LU decomposition that Meta.Numerics also supports. The naive approach of setting AI = A.Inverse()
and multiplying AI * b will work, but a numerical analyst will tell you it has poor performance and stability. You could also use QR decomposition, which is very stable but a bit slower than LU decomposition.



Hi ichbin.
First of all thanks for your response. I've been debugging all the "presolving" code, that's why I took so long to report back.
I followed your example with the LU decomposition, but unfortunatelly Visual Studio gives me a "cannot divide by zero" error back :\
In alternative I tried QR decomposition (which is more stable) and the fact is that I actually managed to run the project. Although this one's output is a series of NaN (not a number), which I'm still wondering why... I've been debugging this project
with an example I had already run on Fortran, and even SGESV (LAPACK single precision) can solve this system, with the very same matrix. Any thoughts ?
Thanks in advance!


Coordinator
Nov 30, 2012 at 8:30 AM

When you say that following my example gives you a DivideByZeroException, I assume you mean this occurs with your matrix, not with the one in my example. The one in my example should work fine, giving you the result x = (1 2 3). First, please verify
and confirm that this works for you.
Next, I will need to see the matrix you are trying to decompose. If it's small enough, just cut and paste it into your reply. If it's too big for that, just let me know; I can send you my email address and you can send it via an attachment in any format
you like.
DivideByZeroException is the expected behavior for a singular matrix, so given only that clue I would just say your problem is not solvable. But if SGESV is solving it something else must be going on. SGESV uses exactly the same algorithm as the Meta.Numerics
example I gave you: LU decomposition with partial pivoting. Possible explanations that come to mind are: (1) there was a problem translating the code to generate your matrix from FORTRAN to VB, so we are not really operating on the same matrix; (2) there is
a bug in LUDecomposition that hasn't shown up on any of the other matrices that have been thrown at it; (3) your matrix is so very close to singular that differences in roundoff between different compilations of the same algorithm have different results.
I look forward to hearing back from you with the matrix, not only to help solve your specific problem, but also to track down any possible bug in Meta.Numerics.


Nov 30, 2012 at 3:34 PM
Edited Nov 30, 2012 at 5:49 PM

Hello again.
Before thinking of sending you the matrix, I decided to try writing the Fortran example matrix into a file which I would load with Visual Studio in order to try to solving it through Meta.Numerics LU decomposition. Altough, once the matrix is really HUGE
(2679 x 2679), the .txt file I made with every element's value is around 280 MB! No way I can load this into Visual Studio without getting a crash. Can the problem be related to the actual size of the matrix? I mean, is there any change I'm running out of
memory somewhere along the calculations ?
Edit: nvm my naiveness. Once I'm debugging with the first iteration results, most of the elements are zero, thus don't need to be declared. I'm doing some testing right now, I'll let you know of the results ASAP. Thanks ! :)



Ok, it seems that I was making two errors.
This A matrix is a stiffness matrix of a structure. From the Fortran program I switched two coordinates on both matrix and vector, which I thought would be innofensive, actually it is not.
The second error was related to the (0) index of the matrix/vector. Once I copied some of the code from the Fortran program, I decided to neglect the (0) index, beggining in (1) as in Fortran. That obviously results in error because this neglection fills
one row and one column of the matrix with zeros.
The LU decomposition is now fully working :) The only problem is that it takes about 1 min to solve each iteration (each LU solve)... On my mac (with a much earlier processor) the Fortran program takes about 2 or 3 seconds per iteration when solving the
very same problem... Is there any chance of speeding this up ?
Thanks again!


Coordinator
Dec 2, 2012 at 7:29 AM

Great, that explains the previous DivideByZeroException, since any matrix with a full row or column of zeros is singular.
Performancewise, can you please clarify what changes in each "iteration"? If only the righthandside vector b changes and the matrix A does not change, you don't need to redo the LU decomposition. In code...
' Don't do this
For Each b As ColumnVector In rightHandSides
Dim LU = A.LUDecomposition()
Dim x = LU.Solve(b)
Next
' Do this instead
Dim LU = A.LUDecomposition()
For Each b As ColumnVector In rightHandSides
Dim x = LU.Solve(b)
Next
I tested this on a dimension2500 system and on my commodity laptop the Solve step took just 200ms. It's the LUDecomposition step that's expensive.
If, on the other hand, your matrix A does change, you will have to redo the LU decomposition. I am pretty skeptical that any (nonsupercomputer) machine could LU decompose such a large matrix in just 23 seconds, even with the Intel MKL installed. That
would be several tens of gigaflops, which is well over the capacity of a typical commodity CPU even if it were doing absolutely nothing else (i.e. no operating system was running) and no memory had to be transfered (i.e. it would have to have a 50MB
L1 cache  a few KB is typical), so I am pretty confident in that skepticism. If however, your system had some special property that allows it to be solved with a less general method than LU decomposition, it's possible that method could be much faster. For
example, if your system is sparse (say 10% or less of the entries in your matrix are nonzero), you might want to try putting it into Meta.Numerics's SparseSquareMatrix and using the iterative Solve method built into that class. For wellbehaved sparse matrices,
an iterative approach can be significantly faster than LU decomposition.

