Linear system solver in D?

[Bill Baxter]

BCS wrote:
> Bill Baxter wrote:
>> Christopher Wright wrote:
>>
>>> BCS wrote:
>>>
>>>> I am going to have a system of equations like this
>>>>
>>>> a_11*x_1 + a_12*x2 + ... a_1n*x_n = y_1
>>>> a_21*x_1 + a_22*x2 + ... a_2n*x_n = y_2
>>>> ..
>>>> ..
>>>> ..
>>>> a_m1*x_1 + a_m2*x2 + ... a_mn*x_n = y_m
>>>>
>>>> y_* and a_* known, I need to find x_*
>>>>
>>>> What is the best available solver for such a system that works under D?
>>>>
>>>> C bindings would work, D code would be better and I'd rather stay 
>>>> away from non portable (uses __ asm and has no port to other system).
>>>>
>>>> If no one knows of a good lib that is ready to use, what is a good C 
>>>> lib to do bindings for?
>>>>
>>>>
>>>> p.s. I'm going to be putting this in a non-linear root finder, has 
>>>> someone already written on of those for D.
>>>
>>>
>>> You definitely want bindings rather than native D. D just hasn't been 
>>> around long enough for people to make decent math libraries for it; 
>>> most of the people with the required skills are still transitioning 
>>> from Fortran.
>>>
>>> You could use GLPK -- it's a linear solver that accepts a superset of 
>>> AMPL. If you're doing serious work on large data sets, go with CPLEX. 
>>> If you manage to write something that does any better than GLPK, 
>>> start a company. CPLEX is significantly better, but you might be able 
>>> to make some money if you marketed it toward smaller research 
>>> projects for $500 or so.
>>
>>
>> Multiarray has bindings to LAPACK.
>>
>> http://www.dsource.org/projects/multiarray
>>
>> There are bindings for GSL which I think uses LAPACK also somewhere on 
>> dsource (either its own project or maybe it was in the 'bindings' 
>> project).
>>
>> I'm working on a new d library that will wrap LAPACK and some sparse 
>> lib like SuperLU and/or TAUCS.  The new lib is based loosely on FLENS.
>>
>> --bb
> 
> thanks, both or you, I'll look at those. Performance isn't that big an 
> issue as I'm only looking at about 15-30 equations and a few minutes run 
> time would be ok, but I'm going to have to run ~1500 passes through it.

Are the equations the same each time?  I.e. do the a_m1 coeffs stay the 
same for all 1500 passes? or do they change each time?

If they stay the same then you can factorize A once and do the much 
faster back-substitution 1500 times.

If you don't care about performance just search the web for some 
C/C++/C#/Java code for Gaussian Elimination with Partial (or better, 
Full) pivoting.  It should be pretty straightforward to port.

Then you won't have to worry about dependencies and building BLAS/LAPACK 
etc.


Actually you didn't say, but are m and n not the same?  If not then you 
need to use least-squares.  If you have more equations than unknowns 
then you cannot generally find an exact solution, but you can find the 
solution that minimizes the 2-norm of the residual ||Ax-b||.  If you 
have too few equations then there are many exact solutions, so generally 
you try to find one that has the smallest 2-norm.  I.e. ||x|| is 
smallest among all possible solutions.

One way to solve a least squares problem is to use the normal equations 
-- you multiply both sides of the equation by A transpose (A^T):
    A^T A x = A^T b

(A^T A) is square, so you can use a regular linear solver on it (like 
gaussian elimination).

--bb