Interesting.As a side note, I thought Julia did support OOP; i.e. multiple dispatch is a form of OOP. I believe Julia's treatment of OOP is similar to CLOS?
maybe, but the absence of member functions in Julia's type makes it different in usage compared to how objects are used in C++ or Python. I don't have a detailed knowledge the technical definitions related to OOP and so I could well be wrong here. Can anyone clarify?
The sparse matrix implementation in julia is fairly young, but it is nice to see that it is already usable. Much of the sparse matrix implementation (except the solvers) is written in julia, and I expect it to get faster as the compiler improves.
In julia, you are using \ to solve the system. This currently defaults to using LU factorization. However, if your system is SPD, there are gains from using a cholesky solver. Matlab does this automatically, and Julia does this automatically in the dense case, but not yet in the sparse case. I suspect that your system is not SPD from your choice of GMRES.
We currently do not have the preconditioned iterative solvers. These should be easy to add, and I suspect this will happen as soon as there is demand.
By default, on 64-bit systems, you get sparse matrices with 64-bit integers. I suspect that there could be a modest performance gain (10-15%) by using sparse matrices with 32-bit integers. Some of this needs to be documented.
I found Amuthan 's blog a while back, but only about two weeks ago I found the time to look seriously at Julia. What I found was very encouraging.
For a variety of teaching and research purposes I maintain a Matlab FEA toolkit called FinEALE. It is about 80,000 lines of code with all the examples and tutorials. In the past week I rewrote the bits and pieces that allow me to run a comparison with Amuthan 's code. Here are the results:
For 1000 x 1000 grid (2 million triangles):
Amuthan's code: 29 seconds
J FinEALE: 86 seconds
FinEALE: 810 seconds
Mind you, we are not doing the same thing in these codes. FinEALE and J FinEALE implement code to solve the heat conduction problem with arbitrarily anisotropic materials. The calculation of the FE space is also not vectorized as in Amuthan's code. The code is written to be legible and general: the same code that calculates the matrices and vectors for a triangle mesh would also work for quadrilaterals, linear and quadratic, both in the pure 2-D and the axially symmetric set up, and tetrahedral and hexahedral elements in 3-D. There is obviously a price to pay for all this generality.
Concerning Amuthan 's comparison with the two compiled FEA codes: it really depends how the problem is set up for those codes. I believe that Fenics has a form compiler which can spit out an optimized code that in this case would be entirely equivalent to the simplified calculation (isotropic material with conductivity equal to 1.0), and linear triangles. I'm not sure about freefem++, but since it has a domain-specific language, it can also presumably optimize its operations. So in my opinion it is rather impressive that Amuthan 's code in Julia can do so well.
Hi Petr: thanks for the note. It is definitely true that making the code more general to accommodate a larger class of elements and/or PDEs would add some computational overhead. I wrote that code primarily to assess the ease with which one could implement a simple FE solver in Julia; the code is admittedly restrictive in its scope and not very extensible either since a lot of the details are hard-coded, but that was a deliberate choice. I had some plans of developing a full fledged FE solver in Julia last year but dropped the idea since I shifted to atomistic simulation and no longer work with finite elements. I would still be interested in hearing about how a general FE solver in Julia performs in comparison to an equivalent FE software. If you have some pointers on that please drop in a note. Thanks!