 Posted Tuesday, August 09, 2005 9:04 PM
|
|
|
|
I have a Matrix A and a Vector B and I want to solve A*X = B. Easy - just do
Matrix X = A.Solve(B);
But I know that B is a vector of zeros with one element equal to one. So I can save some CPU time by doing a LU decomposition. Using Numerical Recipes code this looks like:
Ludcmp(A, indx); // Replaces A by the LU decomposition and sets Vector indx to the row permutation effected by the partial pivoting
Lubksb(A, indx, B); // Solves A*X = B putting the answer in B
This is very efficient as it allows for the case where B has zeros. Can this be done in Bluebit? The Vector indx looks to be the problem
|
|
|
Posted Tuesday, August 09, 2005 11:08 PM
|
|
|
|
Hi,
You can use the LU object of .Net Matrix Library in the same way:
[VB]
Dim luDecA As LU = New LU(A)
[C#]
LU luDecA = new LU(A);
Then you may use the LU decomposition of matrix A with multiple right hand sides stored in matrix B:
X = luDecA.Solve(B)
The permutation matrix can be returned by the P property of the LU object or you can get the pivots in an array by using the the GetPivots method.
You may find more about the LU class here : http://www.bluebit.gr/NET/Library/LU.html
Please let me know if the above has solved your problem.
Trifon Triantafillidis | Lead Developer | 
|
|
|
|
Posted Thursday, August 11, 2005 10:14 PM
|
|
|
|
I know I can use the LU decomposition in the same way. My question was whether or not it would take advantage of the fact that B is a vector of zeros with only one one. Only the Bluebit developer would know that.
My guess is that if it doesn't take advantage of the zeros it would be slower. If it doesn't take advantage of the zeros I could probably write my own back substitution that would.
|
|
|
Posted Friday, August 12, 2005 5:46 AM
|
|
|
|
No it will not make any difference if B is a zero vector.
In order to verify this you may try the following code:
int N = 1000;
Matrix A = new Matrix(N,N);
A.FillRandom();
LU luA = new LU(A);
Matrix B = new Matrix(N,1);
PerformanceTimer pft = new PerformanceTimer();
//Solve with B as a zero vector
pft.StartTimer();
luA.Solve(B);
pft.StopTimer();
Console.WriteLine("Time ={0}", pft.Duration);
// Fill B with random values and solve again
B.FillRandom();
pft.Reset();
pft.StartTimer();
luA.Solve(B);
pft.StopTimer();
Console.WriteLine("Time ={0}", pft.Duration);
In order to run the above example you will need to include in your project the PerformanceTimer class Its code can be found here : http://www.bluebit.gr/forum/topic.asp?TOPIC_ID=58
Trifon Triantafillidis | Lead Developer | |
|
|
|
Posted Wednesday, August 24, 2005 12:09 AM
|
|
|
|
One day I will figure out how to use the GetPivots method to optimise LU code for my special B case. But I don't really have to at the moment. The Bluebit code is an order of magnitude faster than the Numerical Recipes code. So just switching to Bluebit has given me a speedup.
[Note that one reason why the Numerical Recipes code is slower than Bluebit is because the former is written to use C# whereas the latter uses the highly optimised LAPACK code (I presume).]
Thanks, Trifon.
|
|
|
|