David Abrahams wrote:
Yes, it's worth "reinventing" all of those FORTRAN algorithms, because using LAPACK et al. doesn't result in expressive code - every combination of input types has to be handled by a different function with some oddball name. Furthermore, if implementations are not generated programmatically, then they're bug-prone. Even if they are generated programmatically, they are limited (want to work on a matrix of octonions or of small matrices? Sorry -- you're out of luck). Furthermore, pre-generated optimized calculations don't scale. If your FORTRAN library contains a routine to add 3 arrays efficiently, it may not be able to do the same thing for 4 or 5 arrays. With expression templates you can just write the expression and let the metaprograms deal with scaling up the calculation.
I was skimming through Todd Veldhuizen's notes on high performance C++ computing that are accessible from the oonumerics site and I noticed another potential problem with just writing FORTRAN wrappers. According to his notes, FORTRAN compilers are not consistent about symbol naming, so a wrapper needs to know what FORTRAN compiler is being used. The notes are from 2000, so there has been time for the situation to change, but this is an important issue if wrappers are used. John