Thread Subject: Off-topic - how to implement x=a\b in C++?

I am a physiologist and I have used Matlab to develop
computer simulations of working muscles. I am now in the
process of porting the code to C++ so that I can run the
model as a screen-saver on about 100 Windows PCs as a
distributed computing project.

The most time-consuming part of the simulations is
repeatedly solving x=a\b where a is 632 by 632 and >90% of
the non-zero elements are on the tridiagonal lines. I am
currently doing this in C++ using ludcmp() and lubksb()
from Numerical Recipes but the process takes ~25 times
longer than the backslash operator in Matlab.

I have been reading about alternatives to the Numerical
Recipes code and one of the best options seems to be to
call LAPACK routines directly from my C++ code. The trouble
with this is that (probably because I am a physiologist
rather than a computer scientist) it's not obvious how I go
about doing this.

It is not for example yet clear to me how to install LAPACK
on a Windows PC (Visual Studio) or how to call Fortran code
from C++. (One of the advantages of Matlab is that one
doesn't normally have to think about things like this.)

So the bottom line is, can anybody suggest (1) good
tutorials that explain how to install LAPACK with Visual
Studio and call routines from C++ ? or (2) useful
alternatives to Numerical Recipes code for x=a\b with a =
632 x 632 ?

Thanks

Ken

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

Ken,

MATLAB uses many LAPACK in they matrix library. It is quite
seldom that you could beat it easily.

Bruno

"Bruno Luong" <b.luong@fogale.fr> wrote in message <fimic8
$gh1$1@fred.mathworks.com>...
> "Ken Campbell" <campbeks@gmail.com> wrote in message
> <fimhba$2gn$1@fred.mathworks.com>...
> > I am a physiologist and I have used Matlab to develop
> > computer simulations of working muscles. I am now in
the
> > process of porting the code to C++ so that I can run
the
> > model as a screen-saver on about 100 Windows PCs as a
> > distributed computing project.
> >
> > The most time-consuming part of the simulations is
> > repeatedly solving x=a\b where a is 632 by 632 and >90%

<SNIP>

> Ken,
>
> MATLAB uses many LAPACK in they matrix library. It is
quite
> seldom that you could beat it easily.
>
> Bruno

Hi Bruno,

Thanks for the quick reply. I realize that it will be hard
to solve simultaneous equations faster than Matlab but
that's not what I am trying to do.

Instead, I want to obtain a solution to a problem that
could be written in Matlab as x=a\b without using Matlab.

The reason I need to do this is that the final program is
going to be run as a distributed computing project on ~100
PCs. I can't use the Matlab Distributed Computing toolbox
for this because the PCs are in different geographic
locations.

Ken
 

Hi,

do the functions from numerical recipies take care of the
special structure of your matrix? Or do they solve it like
a full matrix? The algorithm is having a complexity of n^3,
so it naturally takes a while.

Do the matrix A in A*x=b change in every iteration?

You could also take a look for the cholesky decomposition.
It applies for positive definite matrices.

I didn't find an good explanation for installing lapack on
windows. Would be interesting to me as well.

REgards,
StefaN


"Ken Campbell" <campbeks@gmail.com> wrote in message
<fiml1j$oqi$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.fr> wrote in message <fimic8
> $gh1$1@fred.mathworks.com>...
> > "Ken Campbell" <campbeks@gmail.com> wrote in message
> > <fimhba$2gn$1@fred.mathworks.com>...
> > > I am a physiologist and I have used Matlab to develop
> > > computer simulations of working muscles. I am now in
> the
> > > process of porting the code to C++ so that I can run
> the
> > > model as a screen-saver on about 100 Windows PCs as a
> > > distributed computing project.
> > >
> > > The most time-consuming part of the simulations is
> > > repeatedly solving x=a\b where a is 632 by 632 and
>90%
>
> <SNIP>
>
> > Ken,
> >
> > MATLAB uses many LAPACK in they matrix library. It is
> quite
> > seldom that you could beat it easily.
> >
> > Bruno
>
> Hi Bruno,
>
> Thanks for the quick reply. I realize that it will be
hard
> to solve simultaneous equations faster than Matlab but
> that's not what I am trying to do.
>
> Instead, I want to obtain a solution to a problem that
> could be written in Matlab as x=a\b without using Matlab.
>
> The reason I need to do this is that the final program is
> going to be run as a distributed computing project on
~100
> PCs. I can't use the Matlab Distributed Computing toolbox
> for this because the PCs are in different geographic
> locations.
>
> Ken
>
>

Hi again,

more that comes to mind. There are blocked LU-
decompositions out there. The complexity goes down
dramatically. The evaluation time as well.
Take a look here:

http://www.osl.iu.edu/research/mtl/tutorial.php3

Is that going to be a boinc project?

Regards,
Stefan


"Stefan " <gonseaw@yahoo.com> wrote in message
<fimmgk$gb9$1@fred.mathworks.com>...
> Hi,
>
> do the functions from numerical recipies take care of the
> special structure of your matrix? Or do they solve it
like
> a full matrix? The algorithm is having a complexity of
n^3,
> so it naturally takes a while.
>
> Do the matrix A in A*x=b change in every iteration?
>
> You could also take a look for the cholesky
decomposition.
> It applies for positive definite matrices.
>
> I didn't find an good explanation for installing lapack
on
> windows. Would be interesting to me as well.
>
> REgards,
> StefaN
>
>
> "Ken Campbell" <campbeks@gmail.com> wrote in message
> <fiml1j$oqi$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.fr> wrote in message
<fimic8
> > $gh1$1@fred.mathworks.com>...
> > > "Ken Campbell" <campbeks@gmail.com> wrote in message
> > > <fimhba$2gn$1@fred.mathworks.com>...
> > > > I am a physiologist and I have used Matlab to
develop
> > > > computer simulations of working muscles. I am now
in
> > the
> > > > process of porting the code to C++ so that I can
run
> > the
> > > > model as a screen-saver on about 100 Windows PCs as
a
> > > > distributed computing project.
> > > >
> > > > The most time-consuming part of the simulations is
> > > > repeatedly solving x=a\b where a is 632 by 632 and
> >90%
> >
> > <SNIP>
> >
> > > Ken,
> > >
> > > MATLAB uses many LAPACK in they matrix library. It is
> > quite
> > > seldom that you could beat it easily.
> > >
> > > Bruno
> >
> > Hi Bruno,
> >
> > Thanks for the quick reply. I realize that it will be
> hard
> > to solve simultaneous equations faster than Matlab but
> > that's not what I am trying to do.
> >
> > Instead, I want to obtain a solution to a problem that
> > could be written in Matlab as x=a\b without using
Matlab.
> >
> > The reason I need to do this is that the final program
is
> > going to be run as a distributed computing project on
> ~100
> > PCs. I can't use the Matlab Distributed Computing
toolbox
> > for this because the PCs are in different geographic
> > locations.
> >
> > Ken
> >
> >
>

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90%
of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The
trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I
go
> about doing this.
>
> It is not for example yet clear to me how to install
LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran
code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?

Not really my area, but I do know that mldivide (\) uses
Gaussian elimination, whereas I _think_ you are attempting
to compute the inverse of your matrix. Others will flame
or add.

Hi,

I think that "ludcmp() and lubksb()" souns like LU-
decomposition, so no matrix inversion.

Regards,
Stefan

"Steve Amphlett" <Firstname.Lastname@Where-I-Work.com>
wrote in message <fimnah$ro4$1@fred.mathworks.com>...
> "Ken Campbell" <campbeks@gmail.com> wrote in message
> <fimhba$2gn$1@fred.mathworks.com>...
> > I am a physiologist and I have used Matlab to develop
> > computer simulations of working muscles. I am now in
the
> > process of porting the code to C++ so that I can run
the
> > model as a screen-saver on about 100 Windows PCs as a
> > distributed computing project.
> >
> > The most time-consuming part of the simulations is
> > repeatedly solving x=a\b where a is 632 by 632 and >90%
> of
> > the non-zero elements are on the tridiagonal lines. I
am
> > currently doing this in C++ using ludcmp() and lubksb()
> > from Numerical Recipes but the process takes ~25 times
> > longer than the backslash operator in Matlab.
> >
> > I have been reading about alternatives to the Numerical
> > Recipes code and one of the best options seems to be to
> > call LAPACK routines directly from my C++ code. The
> trouble
> > with this is that (probably because I am a physiologist
> > rather than a computer scientist) it's not obvious how
I
> go
> > about doing this.
> >
> > It is not for example yet clear to me how to install
> LAPACK
> > on a Windows PC (Visual Studio) or how to call Fortran
> code
> > from C++. (One of the advantages of Matlab is that one
> > doesn't normally have to think about things like this.)
> >
> > So the bottom line is, can anybody suggest (1) good
> > tutorials that explain how to install LAPACK with
Visual
> > Studio and call routines from C++ ? or (2) useful
> > alternatives to Numerical Recipes code for x=a\b with a
=
> > 632 x 632 ?
>
> Not really my area, but I do know that mldivide (\) uses
> Gaussian elimination, whereas I _think_ you are
attempting
> to compute the inverse of your matrix. Others will flame
> or add.

"Stefan " <gonseaw@yahoo.com> wrote in message
<fimmgk$gb9$1@fred.mathworks.com>...
> Hi,
>
> do the functions from numerical recipies take care of the
> special structure of your matrix? Or do they solve it like
> a full matrix? The algorithm is having a complexity of n^3,
> so it naturally takes a while.
>
> Do the matrix A in A*x=b change in every iteration?
>
> You could also take a look for the cholesky decomposition.
> It applies for positive definite matrices.
>
> I didn't find an good explanation for installing lapack on
> windows. Would be interesting to me as well.
>
> REgards,
> StefaN
>

<SNIP>

Hi Stefan,

The Numerical Recipes functions solve the full matrix so
there's probably scope for improvement here.

The (relatively small number of) elements off the
tridiagonal lines in A change after every iteration.

This is very naive but is there an easy way to tell if a
matrix is positive definite? The Wikipedia entry

http://en.wikipedia.org/wiki/Positive_definite_matrix

is tough going for me.

I will let you know if I work how to install LAPACK on Windows.

Ken

"Stefan " <gonseaw@yahoo.com> wrote in message
<fimmgk$gb9$1@fred.mathworks.com>...
> Hi,
>
> do the functions from numerical recipies take care of the
> special structure of your matrix? Or do they solve it like
> a full matrix? The algorithm is having a complexity of n^3,
> so it naturally takes a while.
>
> Do the matrix A in A*x=b change in every iteration?
>
> You could also take a look for the cholesky decomposition.
> It applies for positive definite matrices.
>
> I didn't find an good explanation for installing lapack on
> windows. Would be interesting to me as well.
>
> REgards,
> StefaN
>

<SNIP>

Hi Stefan,

The Numerical Recipes functions solve the full matrix so
there's probably scope for improvement here.

The (relatively small number of) elements off the
tridiagonal lines in A change after every iteration.

This is very naive but is there an easy way to tell if a
matrix is positive definite? The Wikipedia entry

http://en.wikipedia.org/wiki/Positive_definite_matrix

is tough going for me.

I will let you know if I work how to install LAPACK on Windows.

Ken

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>

I do not know well Numerical Recipes. I just want to make a
remark that there is a the there exists quick algorithm for
LU decomposition and resolution for triangular matrix.
Interpolation by 1D cubic spline calls also for such
resolution of triangular matrix, and it is usually pretty fast.

Bruno

"Stefan " <gonseaw@yahoo.com> wrote in message
<fimmgk$gb9$1@fred.mathworks.com>...
> Hi,
>
> do the functions from numerical recipies take care of the
> special structure of your matrix? Or do they solve it like
> a full matrix? The algorithm is having a complexity of n^3,
> so it naturally takes a while.
>
> Do the matrix A in A*x=b change in every iteration?
>
> You could also take a look for the cholesky decomposition.
> It applies for positive definite matrices.
>
> I didn't find an good explanation for installing lapack on
> windows. Would be interesting to me as well.
>
> REgards,
> StefaN
>

<SNIP>

Hi Stefan,

The Numerical Recipes functions solve the full matrix so
there's probably scope for improvement here.

The (relatively small number of) elements off the
tridiagonal lines in A change after every iteration.

This is very naive but is there an easy way to tell if a
matrix is positive definite? The Wikipedia entry

http://en.wikipedia.org/wiki/Positive_definite_matrix

is tough going for me.

I will let you know if I work how to install LAPACK on Windows.

Ken

Ken Campbell wrote:
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?

My first thought here on such a gross time difference and the buzzword
"C/C++" is are you perhaps processing the arrays in "wrong" sequence?
ML stores arrays in column-major order a la its Fortran heritage whereas
C/C++ use row-major storage. If you developed the algorithm under
Matlab and directly moved that to C w/o consideration of this, it
_could_ be the cache miss owing to accessing non-adjacent elements is a
significant problem.

As for LAPACK and C/C++, I've no clue -- I'm a Fortran guy, but for
Fortran it's simply either downloading a platform-specific binary if you
can find one or compiling the routines and linking as an external
library. I wouldn't think there's any real "installation" required;
simply tell VS where the external library is located for the linker.

As Steve says, if I'm way off base somebody will flame/comment... :)

--

"Stefan " <gonseaw@yahoo.com> wrote in message
> Hi again,
>
> more that comes to mind. There are blocked LU-
> decompositions out there. The complexity goes down
> dramatically. The evaluation time as well.
> Take a look here:
>
> http://www.osl.iu.edu/research/mtl/tutorial.php3
>
> Is that going to be a boinc project?
>
> Regards,
> Stefan
>

Hi Stefan,

Sorry for the multiple posts above.

No, it's not a boinc project - I thought about it but ended
up creating my own distributed computing system. It was a
lot of fun and not too hard in the end.

I'm now trying to follow up on your linear algebra tips...

Ken

You could check for positive definiteness with just
applying the chol command inside matlab to your matrix.
If it works for one matrix, you need to check it works for
all matrices.
In my first post I was kind a quick, I would suggest you to
stay with LU, because of numerical stability.

Regards,
Stefan


"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimodm$h94$1@fred.mathworks.com>...
> "Stefan " <gonseaw@yahoo.com> wrote in message
> > Hi again,
> >
> > more that comes to mind. There are blocked LU-
> > decompositions out there. The complexity goes down
> > dramatically. The evaluation time as well.
> > Take a look here:
> >
> > http://www.osl.iu.edu/research/mtl/tutorial.php3
> >
> > Is that going to be a boinc project?
> >
> > Regards,
> > Stefan
> >
>
> Hi Stefan,
>
> Sorry for the multiple posts above.
>
> No, it's not a boinc project - I thought about it but
ended
> up creating my own distributed computing system. It was a
> lot of fun and not too hard in the end.
>
> I'm now trying to follow up on your linear algebra tips...
>
> Ken

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken


Is it the least square problem you mean?

Ax=b

multiply with A' (transpose) on both sides and solve the nxn
linear system:

(A'*A)*x=(A'*b)
x=(A'*A)/(A'*b)% similar to (but faster) inv(A'*A)*(A'*b)

this is the same as x=A\b (backslash)

/Per

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimodm$h94$1@fred.mathworks.com>...

>
> I'm now trying to follow up on your linear algebra tips...
>

Ken,

I strongly recommend you to look at algorithm for
tri-diagonal matrix. The resolution of such system comprises
just two for-loops (one forward, other backward) on the
dimension of your matrix. It is very fast.

Bruno

"Ken Campbell" <campbeks@gmail.com> wrote in message
>
> This is very naive but is there an easy way to tell if a
> matrix is positive definite?
>

Yes, Cholesky as Stephan said. But before this step, please
make sure your matrix is symmetric (upper and lower
diagonals are equal). Is it the case?

Bruno

"Per Sundqvist" <per.sundqvist@uam.es> wrote in message
<fimrkb$1jq$1@fred.mathworks.com>...

<SNIP>
 
 
> Is it the least square problem you mean?
>
> Ax=b
>
> multiply with A' (transpose) on both sides and solve the nxn
> linear system:
>
> (A'*A)*x=(A'*b)
> x=(A'*A)/(A'*b)% similar to (but faster) inv(A'*A)*(A'*b)
>
> this is the same as x=A\b (backslash)
>
> /Per

Hi Per,

I'm not probably not good enough at algebra to understand
this point.

Is x=(A'*A)/(A'*b) a general (fast) way of solving for x?

I also think that I don't have a least-squares problem
because I have the same number of equations as unknowns.

Ken

Hi,

no you don't have a least squares problem.

x = inv(A'*A)*(A'*b)

where inv(A'*A)*A' is called the pseudo inverse.
It includes a matrix-matrix-multiplication, an inverse and
two matrix-vector-multiplications.

I really don't think this one is faster!

Stay with LU.

Regards,
Stefan

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimsri$h5j$1@fred.mathworks.com>...
> "Per Sundqvist" <per.sundqvist@uam.es> wrote in message
> <fimrkb$1jq$1@fred.mathworks.com>...
>
> <SNIP>
>
>
> > Is it the least square problem you mean?
> >
> > Ax=b
> >
> > multiply with A' (transpose) on both sides and solve
the nxn
> > linear system:
> >
> > (A'*A)*x=(A'*b)
> > x=(A'*A)/(A'*b)% similar to (but faster)
inv(A'*A)*(A'*b)
> >
> > this is the same as x=A\b (backslash)
> >
> > /Per
>
> Hi Per,
>
> I'm not probably not good enough at algebra to understand
> this point.
>
> Is x=(A'*A)/(A'*b) a general (fast) way of solving for x?
>
> I also think that I don't have a least-squares problem
> because I have the same number of equations as unknowns.
>
> Ken
>

"Ken Campbell" <campbeks@gmail.com> wrote in message
news:fimhba$2gn$1@fred.mathworks.com...

*snip*

> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?

I only have an answer for the first part of your first question: check out
question 1.13 in the LAPACK FAQ, linked from the main LAPACK page:

http://www.netlib.org/lapack/

Question 1.22 might help you find resources that tell you how to call LAPACK
from C++.

--
Steve Lord
slord@mathworks.com

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

Sounds pretty sparse to me. I'm the author of x=A\b when A
is sparse and square, and you can find all my codes (in C,
and callable from C++) on my web page or on the File Exchange.

Drop me an email with your matrix (davis@cise.ufl.edu) and
I'll take a look and give you a recommendation on which of
my solvers to use (CHOLMOD, UMFPACK, KLU, CSparse, LDL, ...
I think that's all of them...).

Also send me a sample right-hand side, please.

Avoid the Numerical Recipes ... they are awful when it comes
to accuracy considerations and efficiency.

Tim Davis

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

Sounds pretty sparse to me. I'm the author of x=A\b when A
is sparse and square, and you can find all my codes (in C,
and callable from C++) on my web page or on the File Exchange.

Drop me an email with your matrix (davis@cise.ufl.edu) and
I'll take a look and give you a recommendation on which of
my solvers to use (CHOLMOD, UMFPACK, KLU, CSparse, LDL, ...
I think that's all of them...).

Also send me a sample right-hand side, please.

Avoid the Numerical Recipes ... they are awful when it comes
to accuracy considerations and efficiency.

Tim Davis

"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

Sounds pretty sparse to me. I'm the author of x=A\b when A
is sparse and square, and you can find all my codes (in C,
and callable from C++) on my web page or on the File Exchange.

Drop me an email with your matrix (davis@cise.ufl.edu) and
I'll take a look and give you a recommendation on which of
my solvers to use (CHOLMOD, UMFPACK, KLU, CSparse, LDL, ...
I think that's all of them...).

Also send me a sample right-hand side, please.

Avoid the Numerical Recipes ... they are awful when it comes
to accuracy considerations and efficiency.

Tim Davis

"Ken Campbell" <campbeks@gmail.com> wrote in message
...
> Thanks for the quick reply. I realize that it will be hard
> to solve simultaneous equations faster than Matlab but
> that's not what I am trying to do.

It's quite possible to beat MATLAB ... or at least tie it,
since much of the x=A\b code is open source. You can beat
MATLAB backslash because backslash has to figure out what
the matrix properties are and then choose LU, CHOL, QR, etc.
 If you know already, you can skip that can call the right
code directly. Like linsolve in MATLAB ... except linsolve
doesn't work for sparse matrices.

dpb <none@non.net> wrote in message <

...

> My first thought here on such a gross time difference and
the buzzword
> "C/C++" is are you perhaps processing the arrays in
"wrong" sequence?
> ML stores arrays in column-major order a la its Fortran
heritage whereas
> C/C++ use row-major storage. If you developed the
algorithm under
> Matlab and directly moved that to C w/o consideration of
this, it
> _could_ be the cache miss owing to accessing non-adjacent
elements is a
> significant problem.
>
> As for LAPACK and C/C++, I've no clue -- I'm a Fortran
guy, but for
> Fortran it's simply either downloading a platform-specific
binary if you
> can find one or compiling the routines and linking as an
external
> library. I wouldn't think there's any real "installation"
required;
> simply tell VS where the external library is located for
the linker.
>
> As Steve says, if I'm way off base somebody will
flame/comment... :)
>
> --

begin flame
WHOOOOOOSH!
end flame ;-)

Seriously ... the reason his C/C++ code is slow is probably
because he's using a dense solver (in C) in Numerical
Recipes. It's not because of row-major vs col-major. It's
because of dense vs sparse.

Numerical Recipes doesn't have a good sparse solver (it
doesn't have a good dense solver either, for that matter...).

I write in C, and x=A\b is mostly in C. I store arrays by
column, in C, just like they "should" be :-) being an old
Fortran hack myself.

Sparse matrices in MATLAB are stored by column as well ...
even though most of the codes that operate on them are in C/C++.

On Nov 29, 8:11 am, "Ken Campbell" <campb...@gmail.com> wrote:
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90% of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I go
> about doing this.
>
> It is not for example yet clear to me how to install LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

Hi Ken,

It seems that the model will serve as a screensaver. Have you thought
about precomputing the sequence of images that the screensaver will
show and make them into a movie? Then the screensaver sinply shows the
movie. This has the additional advantage that if there are processes
running in the background, the screensaver will not drain the
computational resources from those background processes.

Just my $0.02

Jomar

I can give some input on using LAPACK with C. Newer
versions of Matlab ship with the LAPACK library so
hopefully you might have it already, try

matlabroot\extern\lib\win32\microsoft\libmwlapack.lib

for 32 bit machines and MS Visual Studio. There's no need
to install anything, and MATLAB has some help if you search
for LAPACK.

Grab a Lapack manual and choose the Lapack function you
wish to use (probably dependant on the outcome of this
discussion). The book "LAPACK users' guide" is probably
your best bet. I think the function DGETRS is used for A*X
= B problems with general (real double).

Next create a header file for the function you want to
use. Here's an extract of one I created to use DGEEV
(performs an eigen decomposition)

...
#if !defined(_WIN32)
#define dgeev dgeev_
#endif

extern void dgeev(char* LeftEig,
                  char* RightEig,
                  ...,
                  int* Info);

Because of LAPACK's Fortran roots all inputs to all
functions are pointers. So now make sure the LAPACK
library is attached to your project and you should be good
to go!

Cheers,
Andrew






"Ken Campbell" <campbeks@gmail.com> wrote in message
<fimhba$2gn$1@fred.mathworks.com>...
> I am a physiologist and I have used Matlab to develop
> computer simulations of working muscles. I am now in the
> process of porting the code to C++ so that I can run the
> model as a screen-saver on about 100 Windows PCs as a
> distributed computing project.
>
> The most time-consuming part of the simulations is
> repeatedly solving x=a\b where a is 632 by 632 and >90%
of
> the non-zero elements are on the tridiagonal lines. I am
> currently doing this in C++ using ludcmp() and lubksb()
> from Numerical Recipes but the process takes ~25 times
> longer than the backslash operator in Matlab.
>
> I have been reading about alternatives to the Numerical
> Recipes code and one of the best options seems to be to
> call LAPACK routines directly from my C++ code. The
trouble
> with this is that (probably because I am a physiologist
> rather than a computer scientist) it's not obvious how I
go
> about doing this.
>
> It is not for example yet clear to me how to install
LAPACK
> on a Windows PC (Visual Studio) or how to call Fortran
code
> from C++. (One of the advantages of Matlab is that one
> doesn't normally have to think about things like this.)
>
> So the bottom line is, can anybody suggest (1) good
> tutorials that explain how to install LAPACK with Visual
> Studio and call routines from C++ ? or (2) useful
> alternatives to Numerical Recipes code for x=a\b with a =
> 632 x 632 ?
>
> Thanks
>
> Ken

You can download a trial version of Intel's LAPACK
implementation called Math Kernel Library and call DGELS.

There are Fortran,Java and C examples how to do it.

Free version LAPACK is called ATLAS (the interface is the
same but it is more fidly).

I dont think it is worth to try writing LA methods from scrach.

Cheers,

Laszlo

"Andrew " <awbsmith@itee.uq.edu.au> wrote in message
<fiokkd$nso$1@fred.mathworks.com>...
> I can give some input on using LAPACK with C. Newer
> versions of Matlab ship with the LAPACK library so
> hopefully you might have it already, try
>
> matlabroot\extern\lib\win32\microsoft\libmwlapack.lib
>
> for 32 bit machines and MS Visual Studio. There's no need
> to install anything, and MATLAB has some help if you search
> for LAPACK.
>
> Grab a Lapack manual and choose the Lapack function you
> wish to use (probably dependant on the outcome of this
> discussion). The book "LAPACK users' guide" is probably
> your best bet. I think the function DGETRS is used for A*X
> = B problems with general (real double).
>
> Next create a header file for the function you want to
> use. Here's an extract of one I created to use DGEEV
> (performs an eigen decomposition)
>
> ...
> #if !defined(_WIN32)
> #define dgeev dgeev_
> #endif
>
> extern void dgeev(char* LeftEig,
> char* RightEig,
> ...,
> int* Info);
>
> Because of LAPACK's Fortran roots all inputs to all
> functions are pointers. So now make sure the LAPACK
> library is attached to your project and you should be good
> to go!
>
> Cheers,
> Andrew
>
>
>
>
>
>
> "Ken Campbell" <campbeks@gmail.com> wrote in message
> <fimhba$2gn$1@fred.mathworks.com>...
> > I am a physiologist and I have used Matlab to develop
> > computer simulations of working muscles. I am now in the
> > process of porting the code to C++ so that I can run the
> > model as a screen-saver on about 100 Windows PCs as a
> > distributed computing project.
> >
> > The most time-consuming part of the simulations is
> > repeatedly solving x=a\b where a is 632 by 632 and >90%
> of
> > the non-zero elements are on the tridiagonal lines. I am
> > currently doing this in C++ using ludcmp() and lubksb()
> > from Numerical Recipes but the process takes ~25 times
> > longer than the backslash operator in Matlab.
> >
> > I have been reading about alternatives to the Numerical
> > Recipes code and one of the best options seems to be to
> > call LAPACK routines directly from my C++ code. The
> trouble
> > with this is that (probably because I am a physiologist
> > rather than a computer scientist) it's not obvious how I
> go
> > about doing this.
> >
> > It is not for example yet clear to me how to install
> LAPACK
> > on a Windows PC (Visual Studio) or how to call Fortran
> code
> > from C++. (One of the advantages of Matlab is that one
> > doesn't normally have to think about things like this.)
> >
> > So the bottom line is, can anybody suggest (1) good
> > tutorials that explain how to install LAPACK with Visual
> > Studio and call routines from C++ ? or (2) useful
> > alternatives to Numerical Recipes code for x=a\b with a =
> > 632 x 632 ?
> >
> > Thanks
> >
> > Ken
>

<SNIP>

Thank you for all of the responses so far. There have been
some great suggestions and I am trying to work through some
of the ideas to see if I can get them to work.

I've been lurking on this newsgroup for a few years and I
have learned a lot about Matlab and scientific programming
in the process. It's a great resource.

Thanks again.

Ken

LAPACK is far superior to Numerical Recipes, but LAPACK is
not the right method for this matrix. The matrix tril(A)
has 939 nonzeros, dimension 632-by-632, and its Cholesky
factor using the AMD ordering has just 1385 nonzeros (just
446 entries are fill-in).

It's symmetric positive definite.

tic ; x1 = cholmod2 (A,b) ; toc
tic ; x2 = cs_cholsol (A,b) ; toc
tic ; x3 = ldlsparse (A,[ ],b) ; toc
tic ; x4 = A\b ; toc
tic ; x5 = full(A)\b ; toc

gives elapsed times of:

Elapsed time is 0.001304 seconds.
Elapsed time is 0.001005 seconds.
Elapsed time is 0.000725 seconds.
Elapsed time is 0.001112 seconds.
Elapsed time is 0.063213 seconds.

on my Thinkpad T41 laptop (Pentium M). For the curious, the
matrix is at http://www.cise.ufl.edu/~davis/Campbell.

I'm not sure why x=A\b is a teeny bit faster than the
cholmod2 mexFunction, since x=A\b is calling CHOLMOD in this
case. This is with MATLAB 7.5. cs_cholsol and ldlsparse
are both faster than x=A\b, just because they're lean and
don't need to decide what algorithm to use, as x=A\b does.

The last line, x=full(A)\b, uses LAPACK (plus a teeny amount
of time for full(A)).

These 3 functions are all in the SuiteSparse on the File
Exchange. If you'd like a simple package with lots of
features, use CSparse (CXSparse for 64bit machines).

Oops. The link is

http://www.cise.ufl.edu/~davis/Campbell.mat

sorry ...