## digitalmars.D.learn - Linear system solver in D?

• BCS (15/15) Feb 18 2008 I am going to have a system of equations like this
• Christopher Wright (10/32) Feb 18 2008 You definitely want bindings rather than native D. D just hasn't been
• Bill Baxter (8/43) Feb 18 2008 Multiarray has bindings to LAPACK.
• BCS (4/56) Feb 19 2008 thanks, both or you, I'll look at those. Performance isn't that big an
• Bill Baxter (23/81) Feb 19 2008 Are the equations the same each time? I.e. do the a_m1 coeffs stay the
• BCS (9/33) Feb 19 2008 no such luck. This will be part of a newton Newton–Raphson non-linear ...
• Guillaume B. (4/62) Feb 19 2008 lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it ca...
• Bill Baxter (14/75) Feb 19 2008 Hmm, I was assuming he was after D code. If not, for instance if you
• Guillaume B. (4/83) Feb 19 2008 Well, lp_solve is in C so the bindings should be fairly easy to write.....
• BCS (4/9) Feb 19 2008 how is it's API? What I'll have is an array of real's in memeory.
• Bill Baxter (12/26) Feb 19 2008 Guiallumes link looks to be for integer programming, not floating point
• Guillaume B. (23/53) Feb 20 2008 It actually supports floating points number... And if your C compiler
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.
Feb 18 2008
Christopher Wright <dhasenan gmail.com> writes:
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.
Feb 18 2008
Bill Baxter <dnewsgroup billbaxter.com> writes:
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
Feb 18 2008
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.
Feb 19 2008
Bill Baxter <dnewsgroup billbaxter.com> writes:
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
Feb 19 2008

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

Feb 19 2008
"Guillaume B." <guillaume.b.spam sympatico.ca> writes:
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
Feb 19 2008
Bill Baxter <dnewsgroup billbaxter.com> writes:
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
Feb 19 2008
"Guillaume B." <guillaume.b.spam sympatico.ca> writes:
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
Feb 19 2008

lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it

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?
Feb 19 2008
Bill Baxter <dnewsgroup billbaxter.com> writes:
BCS wrote:

lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it

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?
Guiallumes link looks to be for integer programming, not floating point calcs. LAPACK or something that wraps it is what you want. I highly recommend you look at the BLASLAPACK dir of dsource/projects/multiarray. The LAPACK wrappers have been tested on both Windows and Ubuntu. I have instructions in both dirs for how to install in either case, and there's a simple test program there that does a linear solve on a small matrix. http://www.dsource.org/projects/multiarray/browser/trunk/multiarray/BLASLAPACK/lapacktest.d The blas and lapack are separate from the rest of multiarray, so you can just check out that BLASLAPACK dir and be good to go. --bb
Feb 19 2008
"Guillaume B." <guillaume.b.spam sympatico.ca> writes:
Bill Baxter wrote:

BCS wrote:

lp_solve is nice too: http://lpsolve.sourceforge.net/5.5/ ... and it