|
|
Blitz Support : |
From: Fernando Perez (fperez_at_[hidden])
Date: 2003-10-17 11:40:11
Jorn Attermann wrote:
> I perform extensive numerical simulations. Until now I have used
> higher-level interpreted languages like R and Ox, but for the project I am
> currently working on, these are much too slow (or my implementation is too
> ambitious!). Since my project is very well suited to OOP I am looking for a
> C++ linear algebra package which has vectors, matrices and 3-dimensional
> matrices and the ability to convert between these. Of course, I also need LU
> and QR and other standard linear algebra methods. But above all, the package
> should allow to implement efficient code.
>
> I have been using TNT (with Array1D, Array2D and Array3D) for a while but
> unfortunately this package lacks higher-level constructs, and I do not feel
> sufficiently capable of implementing my own methods.
>
> Before I make a change to another package I would like to ask you: Have you
> had experiences with different C++ linear algebra packages and if so: what
> are the pros and cons. Perhaps you can direct me to some information about
> it.
While I can't comment specifically on interfacing Blitz with existing linear
algebra libraries (others here are better qualified than me for that), let me
suggest you can use a hybrid approach.
I am very happy these days doing most of my work in python, and letting blitz
manage the low-level bottlenecks. The nice thing about this is that Blitz
arrays are essentially identical to Python's Numeric arrays, and the cost of
converting between the two is very small. You can do it yourself by writing
extension modules, but in the early stages of prototyping you don't even need
this much. Using weave.inline() (from the scipy library, http://scipy.org)
you can drop C++ code straight into your Python sources, and use C++ as an
'interpreted' language.
Here's a simple example of a python function with C++ code inside it. The
gnod, gwei, f2cmat and c2fmat variables are Numeric arrays, which the C++ code
transparently sees as Blitz arrays of their respective dimensionalities:
def trmats(gnod,gwei):
"""trmats(gnod,gwei) -> f2cmat,c2fmat
Construct the translation matrices to go back and forth between function
values and interpolating coefficients. """
nnod = len(gnod)
f2cmat = N.zeros((nnod,nnod),'d')
c2fmat = N.zeros((nnod,nnod),'d')
code = """
// Construct transform matrices for normalized polynomials.
#define MAX_NOD 128
double p[MAX_NOD];
double x;
// First build f2cmat
p[0] = 1.0;
for (int j=0;j<nnod;++j) {
x = gnod(j);
p[1] = x;
for (int k=2;k<nnod;++k) {
p[k] = ((2.0*k-1.0)*x*p[k-1] - (k-1.0)*p[k-2])/k;
}
for (int i=0;i<nnod;++i) {
f2cmat(i,j) = sqrt(i+0.5)*p[i]*gwei(j);
}
}
// Now build c2fmat
for (int i=0;i<nnod;++i)
for (int j=0;j<nnod;++j)
c2fmat(i,j) = f2cmat(j,i)/gwei(i);
"""
inline(code,['gnod','gwei','nnod','f2cmat','c2fmat'],
type_converters = converters.blitz)
return f2cmat,c2fmat
And in http://windom.colorado.edu/~fperez/python/python-c you can find a bit
more background on these ideas.
As I said, for maximum performance you can write your own extension modules by
hand (weave does have some overhead), but initially this inlined approach for
prototyping is extremely effective and easy on the developer.
This allows you to retain the fluidity of a high-level language, while easily
avoiding the cost of interpretation only where it makes a difference. With
careful use of this mix, you can get performance which is indistinguishable
from a pure C++ code with a fraction of the development effort.
Best,
f.