From: Dag Lindbo (dag_at_[hidden])
Date: 2008-08-11 05:31:19
Paul C. Leopardi wrote:
> On Thu, 7 Aug 2008, Jesse Perla wrote:
>> I am trying to patch together a set of C++ libraries to provide a really
>> straightforward C++ development toolkit /set of libraries for economists to
>> use as an alternative to fortran/matlab. This is for a set of tool
>> proposals that a number of people might follow for separate projects.
>>
>> I am sure this request has been repeated a million times, but I looked
>> through mailing lists, etc. and cannot seem to find definitive wisdom about
>> what is the best approach to take. Can anyone help out with ideas about
>> where the industry is going? As far as I can see, the basic choice is
>> between boost::ublas and blitz++.
>
> Dear Jesse,
> I went through the same exercise for my GluCat library in 2001-2002. Time and
> C++ have moved on since, but essentially, at the time, the choice was between
> Blitz and MTL.
> http://www.osl.iu.edu/research/mtl/
> MTL won because Blitz was not a matrix library but rather an array and tensor
> library and it was not obvious how to do what GluCat needed, which was
> primarily sparse matrix algebra, including inverses. The MTL mailing list
> archive has been taken down, so I'm including the evaluation at the end of
> this message.
>
> After I implemented GluCat using MTL, Joerg Walter told me about uBLAS and it
> had such obvious advantages that I switched GluCat to using uBLAS.
> http://sourceforge.net/forum/forum.php?forum_id=239344
> See also thread http://listarchives.boost.org/Archives/boost/2002/05/29119.php
> and http://osdir.com/ml/lib.mtl.devel/2002-12/msg00013.html
>
> As I said though, time has moved on. There is now an alpha version of MTL 4
> which may be worth a look, if it is still being actively developed:
> http://www.osl.iu.edu/research/mtl/mtl4/
Jesse,
Paul pretty much covers it. I've been working a lot recently with MTL4
and it is definitely a very strong contender. It has a significant
performance edge over other packages and is very clean to work with.
That said, it is still in development and bugs pop up now and then. It
takes a bit of C++ skill to be able to spot what is a bug in the library
as opposed to a mistake by the user. Also, it is research per se and has
seen a lot less use than uBLAS. If your intended audience does not have
much of C++ background I would recommend uBLAS, but keep an eye on MTL4
as it matures.
/Dag Lindbo
>
>
> ---------- Forwarded Message ----------
>
> Subject: Evaluation of MTL for GluCat
> Date: Sun, 2 Dec 2001
> From: "Paul C. Leopardi" <leopardi_at_[hidden]>
> To: mtl-devel_at_[hidden]
>
> Hi all,
> I have decided to write my critique of MTL as an evaluation of MTL in
> relation to the requirements of GluCat. This may make it easier to understand
> the pros and cons. I'm especially interested to know if I have gotten
> anything wrong. I may also post a longer evalution of POOMA, Blitz++ and MTL
> to http://glucat.sourceforge.net
> Thanks
>
> Evaluation of MTL for GluCat, 2001-12-02
> ----------------------------------------
>
> GluCat requirements for linear algebra libraries in decreasing order of
> importance.
> 0) Basics
> It is assumed that the library is C++ and open source.
> 1) Linear algebra
> GluCat needs to calculate inverses and solve linear systems.
> 2) Sparsity
> GluCat needs to use sparse storage for both space and time savings, especially
> for basis generation.
> 3) Accuracy and stability of linear algebra
> Glucat must accurately convert between the coordinate representation and the
> matrix representation.
> 4) Current support
> GluCat will be publically available and will have a public dependence on the
> linear algebra library.
> 5) Simple semantics for assignment to uninitialized matrices
> GluCat needs to use structures containing matrices where the shape is not
> known in advance.
> 6) Simple syntax for linear operations: matrix algebra notation
> GluCat algorithms view both matrix-matrix and matrix-vector multiplication as
> primitive operations. Syntax should be similar to Matlab.
> 7) Simple semantics for linear operations
> GluCat algorithms should be expressed as simply as possible with no
> pre-conditions, special cases or side effects. Semantics
> should be similar to Matlab.
>
> Summary of evaluation of MTL
> PROS
> 1) Linear algebra including LU solvers
> 2) Support for sparsity
> CONS
> 1) Cannot use copy() with uninitialized variables
> Does not meet requirement 5
> 2) No matrix algebra syntax
> Does not meet requirement 6
> 3) Need to pre-allocate storage for sparse matrix operations
> Does not meet requirement 7
> 4) Needs more sophisticated linear algorithms for accuracy and stability
> Does not meet requirement 3
>
> Decision:
> Use MTL for GluCat, with the following workarounds:
>
> 1) Cannot use copy() with uninitialized variables
> Avoid copy() where possible, initialize variables where possible, and use
> B = matrix(A.nrows(), A.ncols()) before copy(A,B) otherwise.
> 2) No matrix algebra syntax
> Use MTL syntax, with Matrix algebra syntax in comments.
> 3) Need to pre-allocate storage for sparse matrix operations
> Use a modified version of MTL which has auto-resize within sparse multiply.
> The alternative would have been to create a temporary with sufficient space
> at every point where a matrix multiplication is done.
> 4) Needs more sophisticated linear algorithms for accuracy and stability
> Add iterative refinement at the point where lu_solve is used.
>
> Details
> PROS
> 1) Linear algebra including LU solvers
> 2) Support for sparsity
>
> CONS
> 1) Cannot use copy() with uninitialized variables.
> MTL splits assignment functionality into a "shallow" operator= and a "deep"
> copy(). Unfortunately, copy() only copies content and not shape. It assumes
> that the matrix being copied to has already been initialized.
>
> Related Bugs:
> 1.1) Documentation should describe this behaviour of copy() but does not
> You either need to expect this as "normal" behaviour, or discover this
> behaviour by using the debugger and examining the source code.
>
> 2) No matrix algebra syntax
> Lack of a matrix algebra syntax makes it necessary to translate matrix
> algorithms into MTL code, which resembles writing a program in assembly
> language. This is time consuming and error prone.
>
> Example
> // r = A*x - b;
> Vector_t r(A.nrows(), 0.0);
> mtl::mult(A, x, r);
> mtl::add(mtl::scaled(b, -1.), r);
>
> 3) Need to pre-allocate storage for sparse matrix operations
> If you want to create sparse matrices you need to know how much storage they
> will need before you operate on them. This is in contrast to eg. Matlab,
> where a sparse array is extended as needed.
>
> Related bugs:
> 3.1) Sparse matrix multiply assigns 0 when it runs out of space
> Sparse multiply outputs a cryptic message and assigns 0 if not
> enough storage has been assigned to sparse array. Either the message should
> be more informative or there should be an option to tell operations to
> automatically resize, such as an auto-resize policy class.
> 3.2) Documentation should mention this behaviour of multiply and does not
> You either need to expect this as "normal" behaviour, or discover this
> behaviour by using the debugger and examining the source code.
>
> 4) Needs more sophisticated linear algorithms for accuracy and stability
> For example, lu_solve uses Gaussian elimination with partial pivoting, but
> does not use iterative refinement.
>
> See also MTL TODO List:
> http://www.osl.iu.edu/research/mtl/todo_list.php3
>
> On Wed, 28 Nov 2001 09:55, Paul C. Leopardi wrote:
>> Hi all,
>> I have been using MTL as the linear algebra engine for the GluCat library (
>> http://glucat.sourceforge.net ) as a project for my Masters by coursework
>> at UNSW. I am currently writing up my project report. One section will be a
>> review and critique of MTL.
>>
>> In the critique of MTL, I will review some of the pitfalls I encountered in
>> starting from a background which is more used to Matlab and Maple and then
>> trying to apply this to writing algorithms using MTL. The main point will
>> be that MTL forces its users to write code at a level lower than Matlab,
>> Fortran 90 or even Blitz++, and that this creates many opportunities for
>> errors and confusion.
>>
>> I intend to publish an early draft on this mailing list, to give other MTL
>> users an opportunity to set me straight.
>> Best regards
>
>
> ---------- Forwarded Message ----------
>
> Subject: Re: MTL critique
> Date: Sun, 3 Feb 2002
> From: "Paul C. Leopardi" <leopardi_at_[hidden]>
> To: szummer_at_[hidden]
>
> On Sun, 3 Feb 2002 10:59, you wrote:
> Hi Martin,
> At this point I don't have much more than my evaluation of MTL versus the
> requirements of GluCat,
> http://www.osl.iu.edu/MailArchives/mtl-devel/msg00244.php
> http://www.osl.iu.edu/MailArchives/mtl-devel/msg00246.php
> but I will attempt to answer your questions below.
> Best regards
>
>> HI Paul,
>>
>> I read your preliminary critique of MTL on the MTL-devel mailing list,
>> written 2001-12-02.
>> Great job! I wish there were more work like yours.
>> You mentioned you will post a more complete critique later - could you
>> CC it to me, please?
>>
>> I would be interested in a critique containing the following (I am happy
>> to wait for your final review, or if you know the answers to some of
>> these questions already, I would be very interested!)
>
>> * performance of MTL compiled with both gcc and KAI - I have only seen
>> the authors' own benchmarks; how about an independent benchmark on the
>> speed of MTL, compared to low-level libraries like ATLAS, or
>> vendor-BLAS?
> I do not have access to the KAI compiler and currently don't have time to
> benchmark against LAPACK and BLAS - MTL is actually at a level in between
> LAPACK and BLAS. While performance matters with GluCat, structure and
> correctness were far more important for the initial build stage.
>
>> * comparison of MTL with Blitz++, TNT, LinAl, LinAlg, and other
>> competing C++ libraries.
> Against Blitz++: Blitz++ does not support general sparse matrices and does
> not provide linear algebra routines such as lu_solve, etc.
> Against TNT: From what I saw of the TNT web site when I examined it a year
> ago, TNT looked unsupported, or at least not currently maintained. I could be
> wrong TNT does not support general sparse matrices.
> I am not familiar with LinAl or LinAlg.
>
>> * practical experience of using MTL - I have heard that it can be
>> difficult to decipher compiler errors from the many nested template
>> expansions
> Deciphering compiler errors is initially difficult but gets much easier very
> quickly. Start with small pieces of code. It is not compiler errors, but
> runtime errors which are the real horror. Compile with debugging on (-g3 for
> g++) and make sure you have a good debugger which can find its way through
> all the long mangled names.
>
>> * how does MTL foster bug-free research code: i) rate of bugs in MTL
>> itself ii) how easy it is to create buggy code in MTL (due to weak
>> types, no array bounds checking etc)
> The low level of the library plus the "shallow copy" semantics means that it
> is easy to create buggy code. Luckily most bugs show up as segmentation
> violations.
>
> You need to be very careful in writing code using MTL. For example, if MTL
> had operators ( such as *= ) and value semantics, I could write:
> matrix_t result = unit(dim);
> const matrix_t* e = Some_Vector_of_Matrices;
> for (index_t k = folded_min; k <= folded_max; ++k)
> if (folded_set[k])
> result *= e[k];
>
> instead of the following:
> result = matrix_t(dim, dim);
> unit(dim, result);
> const matrix_t* e = Some_Vector_of_Matrices;
> for (index_t k = folded_min; k <= folded_max; ++k)
> if (folded_set[k])
> {
> matrix_t temp(dim, dim);
> mtl::mult(result, e[k], temp);
> mtl::copy(temp, result);
> }
>
> With MTL, you need to use mtl::copy to prevent aliasing, but when you do, you
> must make sure that the matrix has already been given a size.
> Because of the ever present likelihood of aliasing, I programmed very
> defensively and probably ended up using mtl::copy more that was strictly
> necessary.
>
>> From a practical point of view: would you recommend using MTL right now?
>> I am trying to decide what library to use.
> It very much depends on your requirements. Mine were:
> 1) Must be a C++ template library
> 2) Must provide linear algebra operations
> 3) Must support general sparse matrices
> MTL fit these requirements and Blitz++, TNT did not.
> POOMA and deal.II were too big to evaluate quickly.
>
>> Thanks,
>>
>> -- Martin
> ---
> _______________________________________________
> ublas mailing list
> ublas_at_[hidden]
> http://listarchives.boost.org/mailman/listinfo.cgi/ublas