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.

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.

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

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.

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

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.


 lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can
read many file formats...

 Guillaume

Guillaume B. wrote:
> 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.
> 
> 
>  lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can
> read many file formats...

Hmm, I was assuming he was after D code.  If not, for instance if you just have some equations you want to solve, then I'd just write a short python/numpy script:
---
from numpy import *
A = asarray([[a_11, a_12 ...],
             [a_21, a_22 ...],
             ...
             [a_m1, a_m2 ...]])
y = asarray([y_1, y_2...])
x = lstsq(A,y)[0]
---

Voila

Bill Baxter wrote:

> Guillaume B. wrote:
>> 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.
>> 
>> 
>>  lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it can
>> read many file formats...
> 
> Hmm, I was assuming he was after D code.  If not, for instance if you just have some equations you want to solve, then I'd just write a short python/numpy script:
> ---
> from numpy import *
> A = asarray([[a_11, a_12 ...],
>               [a_21, a_22 ...],
>               ...
>               [a_m1, a_m2 ...]])
> y = asarray([y_1, y_2...])
> x = lstsq(A,y)[0]
> ---
> 
> Voila

 Well, lp_solve is in C so the bindings should be fairly easy to write... I
suggested that since others suggested using bindings to existing libraries.

 Guillaume

Reply to Guillaume B.,

> lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it
> can read many file formats...
> 
> Guillaume
> 

how is it's API? What I'll have is an array of real's in memeory.

Hmm that's a though, does anyone know if any of the listed libs work with 80bit FP?

Reply to Bill,

> BCS wrote:
> 
>> 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.
> 

no such luck. This will be part of a newton Newton–Raphson non-linear solver for doing a numerical approximation of a system of non-linear differential equations. (a.k.a. bad, worse, nasty)

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

performance takes a back seat to complexity but is still an issue. OTOH generating bindings for canned C code is less complex than porting it.

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

building and what not isn't a big issue.

> Actually you didn't say, but are m and n not the same?

On the way out the door after posting that I relized that I put in NxM rather than NxN. oops.

> 
> --bb
>