@halirutan Hi, halirutan. I use fortran, because I need speed. Before that I have tried my best to write compiled function in mma, but it is 10 times slower. I insist on investigating it because I think mixing fortran is a essential method for people who need extreme speed, and also I saw several post on this site asking for method for calling external numerical fortran library. If this kind of calling is defective, though maybe safe in single thread, it is dangerous.

@halirutan I actually made a post mathematica.stackexchange.com/q/105797/4742

@matheorem OK, but for which algorithms do you need speed? There is absolutely no way that a C program for the small code you provided in Fortran would be any slower.

Contrary, I would claim that due to the latest development in C/C++ compilers, you should get it probably faster.

@halirutan You mean probably faster than fortran? C is definitely not suited for matrix arithmetic. My subroutine is matrix intensive. For c++, there is Eigen library which is handy. But for my subroutine it is still slower than fortran. And the compiler is important, for my fortran subroutine gfortran is much slower than intel fortran, while for my c++ routine, g++ is faster than intel icl. Intel fortran is really an excellent compiler.

@matheorem But since you have the Intel compiler, have you tried to use the C-library like lapack?

@halirutan yeah, but lapack is designed for large matrix, it is slow for small matrix. So I write linearsolve by myself.

@matheorem Exactly, and at which point in this very simple code is Fortran supposed to be faster?

(@matheorem So we are talking about e.g. twobytwolinearsolve, right?)

@halirutan my version in the post is a special treatment of 2x2, using direct formula. According to my googleing, it is a general consensus, for inversion or something else, matrix below 4x4, using direct formula is efficient. For general LU decomposition, a self written version is efficient if the dimension is below 10. here is benchmark I did month ago, postimg.org/image/l9ofl2hkh the green is fortran calling intel lapack, the others are routines I write.

@halirutan I make "LUno-outer" as reference 1

@halirutan the z axes is the time divided by timing of "LUno-outer"

And you are doing the long chain of Mathematica calling C calling Fortran, because you like to write things down in Fortran?

@halirutan Er...I just feel that I maybe need fortran from time to time for speed. Actually I've just started using librarylink. Before that I use only mma for numerically computation, at that time, I calculate Eigenvalues of relatively large matrix, so mathematica is pretty good, since it call MKL.

@halirutan But now the problem change, the reason I can't stand a 10 times slower compiled function of mma is that after I boost 10 times faster using fortran or C++ in the calculation intensive part, I can directly do "just in time" plotting(for this particular problem), don't have to wait too much time, you know that feeling of waiting: ) and even manipulation, that is awesome.

@halirutan About algorithm, I definite not being able to design LU decomposition algorithm, I read books, such as "matrix computation" by Golub, and then write code according the book

@matheorem OK, it's definitely good that you try to implement that yourself and you should go on. Just keep the warning in mind that things like the listable difference can happen and that implementing a numerical scheme is hard. While you can be kind of sure that LAPACK algorithms were tested pretty well, you have to expect your algorithms to fail at some point. Be it because you used a compiler switch that does some weird optimization or be it that you lost some accuracy to rounding errors.

@matheorem And just doing it is the only way of getting deep knowledge at something.

@halirutan yeah, can't agree more