## 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
```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

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

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

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

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
```Reply to Guillaume B.,

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

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

It actually supports floating points number... And if your C compiler
supports 80 bits floating point, you can probably compile it to support it
since the functions are like this:

unsigned char add_constraint(lprec *lp, REAL *row, int constr_type, REAL
rh);

Although by default, REAL is probably a typedef for double.

And lp_solve is actually a solver for Linear Programs... It's an optimizer,
not a simple solver for a system of linear equations... It's to solve
things like this:

maximize 143x + 60y
where
120x + 210y <= 15000
110x + 30y <= 4000
x + y <= 75
x >= 0
y >= 0

So it might not be exactly what you want... You can see an example problem
(the explanations are very good) with an example of the C API there:
http://lpsolve.sourceforge.net/5.5/formulate.htm

Hope this helps!

Guillaume
```
Feb 20 2008