A new Python interface for Elemental

[Michael C. Grant]

Hey folks,

Over the last couple of weeks Jack and I have had a bit of a sprint to finish up a SWIG-based Python interface for Elemental. A preliminary version has been committed to Github. We've already found it to be quite usable, at least for demonstrating and testing purposes.

It's very easy to build and try now, if you have Python with NumPy>=1.7.1 installed. Just add USE_SWIG=ON to your CMake command:

cd elemental
mkdir build
cd build
cmake -D USE_SWIG=ON ..
make

Once it is built, you can run it right from the build directory. (Note: Anaconda/EPD users *may* need to set LD_LIBRARY_PATH or DYLD_LIBRARY_PATH to get this to work.)

python
>>> import elem
>>> A = elem.DistMatrix_d()
>>> elem.Walsh( A, 4 )
>>> elem.Display( A )

I know that SWIG was probably not on the top of anyone's list when we discussed this back in late March. I spent a fair amount of time working with SWIG back in late March and April, but I eventually decided to work on some architectural changes to the C++ code itself (which are not yet part of the main repo). Recently, Jack and I both decided that ultimately something needed to Just Get Done, and that my existing SWIG work was the closest we had to this goal.

I'm afraid we haven't written a lick of documentation for it yet: this email is it. But we haven't modified the calling formats, so elem::FUNC( A, B, C ) in C++ just becomes elem.FUNC( A, B, C ) in Python. Yes, that means it's not very Pythonic at times, such as when the return value is written to one of the arguments, but at least it will make it easy for someone to move between C++ and Python for now.

We welcome your feedback and bug reports! We are well aware of the Python extension is well over twice as large as the C++ Elemental library. It's quite likely there are some serious efficiency problems buried in here, say in some of the more heavily overloaded functions. Please don't hesitate to report those as you encounter them in practice, but note that the solution may involve some of the longer-term work I'm doing on the C++ side. In the short term, I intend to focus on doing a better job of NumPy integration and perhaps exploring more Python-friendly calling sequences.

Regards,
Michael Grant

------

--- You do not need to Initialize or Finalize elemental. Initialization is performed automatically on module load, and finalization on exit. However, if for some reason you really, really *want* to run elem.Finalize() before exit, go for it, it won't bite.

--- The constructors for DistMatrix<T,MC,MR> have a single-letter data type suffix:
elem.DistMatrix_s - float
elem.DistMatrix_d - double
elem.DistMatrix_c - Complex<float>
elem.DistMatrix_z - Complex<double>
Constructors for other DistMatrix<T,U,V> objects append U and V suffixes accordingly; e.g.
elem.DistMatrix_d_VC_STAR - DistMatrix<double,VC,STAR>

--- Functions are overloaded, and do not have distribution or data type suffixes.

--- elem::Matrix objects are wrapped in NumPy arrays. We've tried to be logical about this. For instance,
A = elem.DistMatrix_d()
elem.Walsh( A, 4 )
B = A.Matrix()
B
returns a 2-D NumPy array containing a read-write *view* of the local matrix data. That means you can use standard NumPy operations to operate on this local data, and any changes will be reflected in the DistMatrix.

--- You can still use the non-distributed, elem::Matrix versions of common Elemental operations. For instance, this works:
B = numpy.empty(())
elem.Walsh( B, 3 )
B
However, note that Elemental uses Fortran ordering, and NumPy by default uses C ordering. My code tries to do the correct thing whenever possible, but it's not particularly intelligent about it, so it often comes at a cost of efficiency. For example, consider this code:
B = numpy.ones((8,8))
elem.Scale( 2, B )
B
The first line constructs a matrix with C ordering. When elem.Scale is called, it creates a Fortran-ordered *copy* of B for its own use. When its operation is complete, it copies the results *back* to the original B array, and deletes its temporary. The good news is that if you make sure your NumPy arrays are Fortran-ordered, these shenanigans don't happen, and Elemental can operate directly on the underlying data:
B = numpy.ones((8,8),order='F')
elem.Scale( 2, B )
B
It can even work on NumPy views and slices:
B = numpy.ones((8,8),order='F')
elem.Scale( 2, B[1:5,1:5] )
B

--- TWO CAVEATS for working with NumPy+Elemental:

1) Some NumPy operations, like resizing, are incompatible with views. Those operations will force you to make a copy of the data, which will "disconnect" you from the source DistMatrix. One way to verify that things are working properly is to look at the OWNDATA flag. For instance, consider this:
A = elem.DistMatrix_d()
elem.Walsh( A, 4 )
B = elem.Matrix()
B.flags
If you do this, you'll see OWNDATA : False, because B is a view for the DistMatrix data. If OWNDATA : True, then this is a copy, and it cannot possibly be connected to a DistMatrix. (For NumPy experts, the matrix B here does *not* have a "base" array, so the expression "B.base is None" returns True even if you haven't made a copy yet. I think I need to fix this.)

2) Currently, 2-D NumPy *arrays* are returned, not *matrices*. This means that operations like multiplication are interpreted in the array sense. I intend to change this eventually.

[Jack Poulson]

For those of you that are interested, the talk I gave at ICS on Monday provides at least a bit more documentation on the new Python interface (as well as some recent streamlining of the C++ boiler-plate code):
http://elemental.googlecode.com/files/using-elem.pdf

I hope to release version 0.80 by the end of the week.

Jack

--
You received this message because you are subscribed to the Google Groups "elemental-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to elemental-de...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.