Tag Archives: FEM

Introducing DifferentialEquations.jl

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/introducing-differentialequations-jl/

Differential equations are ubiquitous throughout mathematics and the sciences. In fact, I myself have studied various forms of differential equations stemming from fields including biology, chemistry, economics, and climatology. What was interesting is that, although many different people are using differential equations for many different things, pretty much everyone wants the same thing: to quickly solve differential equations in their various forms, and make some pretty plots to describe what happened.

The goal of DifferentialEquations.jl is to do exactly that: to make it easy solve differential equations with the latest and greatest algorithms, and put out a pretty plot. The core idea behind DifferentialEquations.jl is that, while it is easy to describe a differential equation, they have such diverse behavior that experts have spent over a century compiling different ways to think about and handle differential equations. Most users will want to just brush past all of the talk about which algorithms simply ask: “I have this kind of differential equation. What does the solution look like?”

DifferentialEquations.jl’s User Interface

To answer that question, the user should just have to say what their problem is, tell the computer to solve it, and then tell the computer to plot it. In DifferentialEquations.jl, we use exactly those terms. Let’s look at an Ordinary Differential Equation (ODE): the linear ODE . It is described as the function

y^\prime = \alpha y

To use DifferentialEquations.jl, you first have to tell the computer what kind of problem you have, and what your data is for the problem. Recall the general ordinary differential equation is of the form

y^\prime = f(y)

and initial condition u_0, so in this case, we have an ODE with data f(y)=\alpha y and u_0. DifferentialEquations.jl is designed as a software for a high-level language, Julia. There are many reasons for this choice, but the one large reason is its type system and multiple dispatch. For our example, we tell the machine what type of problem we have by building a DEProblem type. The code looks like this:

using DifferentialEquations
alpha = 0.5 #Setting alpha to 1/2
f(y,t) = alpha*y
u0 = 1.5
prob = ODEProblem(f,u0)

where prob contains everything about our problem. You can then tell the computer to solve it and give you a plot by, well, solve and plot:

timespan = [0,1] # Solve from time = 0 to time = 1
sol = solve(prob,timespan) # Solves the ODE
plot(sol) # Plots the solution using Plots.jl

And that’s the key idea: the user should simply have to tell the program what the problem is, and the program should handle the details. That doesn’t mean that the user won’t have access to to all of the details. For example, we can control the solver and plotters in more detail, using something like

sol = solve(prob,alg=:RK4) # Unrolled and optimzed RK4
plot(sol,lw=3) # All of the Plots.jl attributes are available

However, in many cases a user may approach the problem for which they don’t necessarily know anything about the algorithms involved in approximating the problem, and so obfuscating the API with these names is simply confusing. One place where this occurs is solving stochastic differential equations (SDEs). These have been recently growing in popularity in many of the sciences (especially systems biology) due to their added realism and their necessity when modeling rare and random phenomena. In DifferentialEquations.jl, you can get started by simply knowing that an SDE problem is defined by the functions f and g in the form

dX_t = f(X_t,t)dt + g(X_t,t)dW_t,

with initial condition u_0, and so the steps for defining and solving the linear SDE is

g(u,t) = 0.3u
prob = SDEProblem(f,g,u0)
sol = solve(prob,timespan)
plot(sol)

If you wish to dig into the manual, you will see that the default solver that is used is a Rossler-SRI type of method and will (soon) have adaptivity which is complex enough to be a numerical analysis and scientific computing research project. And you can dig into the manual to find out how to switch to different solvers, but the key idea is that you don’t have to. Everything is simply about defining a problem, and asking for solutions and plots.

And that’s it. For more about the API, take a look at the documentation or the tutorial IJulia notebooks. What I want to discuss is why I believe this is the right choice, where we are, and where we can go with it.

What exactly does that give us?

Julia was created to solve the many-language problem in scientific computing. Before people would have to write out the inner loops as C/Fortran, and bind it to a scripting language that was never designed with performance in mind. Julia has done extremely well as solving this problem via multiple-dispatch. Multiple dispatch is not just about ease of use, but it is also the core of what allows Julia to be fast . From a quote I am stealing from IRC: “Julia: come for the fast loops, stay for the type-system”.

In my view, the many-language problem always had an uglier cousin: the many-API problem. Every package has its own way of interacting with the user, and it becomes a burden to remember how all of them work. However, in Julia there seems to be a beautiful emergence of packages which solve the many-API problem via Julia’s multiple-dispatch and metaprogramming functionalities. Take for example Plots.jl. There are many different plotting packages in Julia. However, through Plots.jl, you can plot onto any “backend” (other plotting package) using just one API. You can mix and match plotting in PyPlot (matplotlib), GR, Plotly, and unicode. It’s all the same commands. Another example of this is JuMP. Its core idea is solver independence: you take your optimization problem, define the model in JuMP’s lingo, and then plug into many different solvers all by flipping a switch.

DifferentialEquations.jl is extending this idea to the realm of differential equations. By using the keyword `alg=:ode45`, the solver can call functions from ODE.jl. And changing it to `alg=:dopri5`, DifferentialEquations.jl will solve your ODE problem using the coveted dopri5 Fortran software. The complexity of learning and understanding many different APIs is no longer a requirement for trying different algorithms.

But why “Differential Equations”? Isn’t that broad?

Sure, there are packages for solving various types of differential equations, all specializing in one little part. But when I was beginning my PhD, quickly found that these packages were missing something. The different types of differential equations that we encounter are not different but many times embody the same problem: a PDE when discretized is a system of ODEs, the probability distribution of evolving SDEs is a PDE (a form of the Heat Equation), and all of the tools that they use to get good performance are the same. Indeed, many contemporary research questions can be boiled down to addressing the question: what happens if we change the type of differential equation? What happens if we add noise to our ODEs which describe population dispersal? What happens if we add to our model that RNA production is reading a delayed signal? Could we make this high-dimensional PDE computationally feasible via a Monte Carlo experiment combined with Feynman-Kac’s theorem?

Yet, our differential equations libraries are separate. Our PDEs are kept separate from our SDEs, while our delay equations hang out in their own world. Mixing and matching solvers requires learning complex APIs, usually large C/Fortran libraries with opaque function names. That is what DifferentialEquations.jl is looking to solve. I am building DifferentialEquations.jl as a hub for differential equations, the general sense of the term.

If you have defined an SDE problem, then via the Forward Kolmorogov equation there is a PDE associated to the SDE. In many cases like the Black-Scholes model, both the SDE and the PDE are canonical ways of looking at the same problem. The solver should translate between them, and the solver should handle both types of differential equations. With one API and the functionality for these contained within the same package, no longer are they separate entities to handle computationally.

Where are we currently?

DifferentialEquations.jl is still very young. Indeed, the project only started a few months ago, and during that time period I was taking 6 courses. However, the package already has a strong core, including

In fact, there are already a lot of features which are unique to DifferentialEquations.jl:

You may have been thinking, “but I am a numerical analyst. How could this program help me?”. DifferentialEquations.jl has a bunch of functionalities for quickly designing and testing algorithms. All of the DEProblems allow for one to give them the analytical solution, and the solvers will then automatically calculate the errors. Thus by using some simple macros, one can define new algorithms in just a few lines of code, test the convergence, benchmark times, and have the algorithm available as an `alg` option in no time (note: all of the ODE solvers were written in one morning!). Thus it is easy to define the loop, and the rest of the functionality will come by default. It’s both a great way to test algorithms, and share algorithms. Contributing will both help you and DifferentialEquations.jl!.

Where are we going?

I have big plans for DifferentialEquations.jl. For example:

  • I will be rolling out an efficiency testing suite so that one could just specify the algorithms you’d like to solve a problem, and along what convergence axis (i.e. choose a few \Delta ts, or along changing tolerances), and it will output comparisons of the computational efficiencies and make some plots. It will be similar in functionality to the ConvergenceSimulation suite.
  • Finite difference methods for Heat and Poisson equation. These are long overdue for the research I do.
  • Changing the tableaus in ODEs and SDEs to StaticArrays so they are stack allocated. This has already been tested and works well on v0.5.
  • Higher-order methods for parabolic SPDEs (a research project with promising results!).
  • Blazing fast adaptivity for SDEs. (Once the paper I have submitted for it is published, it will be available. It’s already implemented!)
  • High-stability high order methods for SDEs (another research project).
  • Parallel methods. I have already implemented parallel (Xeon Phi) solvers and described them in previous blog posts. They simply need to be integrated into DifferentialEquations.jl. I would like to have native GPU solvers as well.
  • Delay and algebraic differential equations.
  • Wrapping more popular solvers. I’d like to add Sundials, LSODE, and PetsC to the list.
  • A web interface via Escher.jl to define DEProblems and get the solution plots. I am looking to have this hosted as an XSEDE Gateway.

If you’d like to influence where this project is going, please file an issue on the Github repository. I am always open for suggestions.

I hope this gives you a good idea on what my plans are for DifferentialEquations.jl. Check out the documentation and give the package a whirl!

The post Introducing DifferentialEquations.jl appeared first on Stochastic Lifestyle.

Julia iFEM3: Solving the Poisson Equation

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem3/

This is the third part in the series for building a finite element method solver in Julia. Last time we used our mesh generation tools to assemble the stiffness matrix. The details for what will be outlined here can be found in this document. I do not want to dwell too much on the actual code details since they are quite nicely spelled out there, so instead I will focus on the porting of the code. The full code will be available soon on a github repository, and so since most of it follows a straight translation from the linked document, I’ll leave it out of the post for you to find on the github.

The Ups, Downs, and Remedies to Math in Julia

At this point I have been coding in Julia for over a week and have been loving it. I come into each new function knowing that if I just change the array dereferencing from () to [] and wrap vec() calls around vectors being used as indexes (and maybe int()), I am about 95% done with porting a function. Then I usually play with cosmetic details. There are a few little details which make code in Julia a lot prettier. For example, for the FEM solver, we need to specify a function. In MATLAB, specifying such the function for which we want to solve -Delta u = f in-line would be done via anonymous functions like:

f = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2));

I tend to have two problems with that code. For one, it’s not math, it’s programming, and so glancing at my equations and glancing at my code has a little bit of a translation step where errors can happen. Secondly, in many cases anonymous functions incur a huge performance decrease. This fact and the fact that MATLAB’s metaprogramming is restricted to string manipulation and eval destroyed a project I had tried a few years ago (general reaction-diffusion solver with a GUI, the GUI took in the reaction equations, but using anonymous functions killed the performance to where it was useless). However, in Julia functions can be defined inline and be first class, and have a nice appearance. For example, the same function in Julia is:

f(x) = sin(2π.*x[:,1]).*cos(2π.*x[:,2]);

Two major improvements here. For one, since the interpreter knows that variables cannot start with a number, it interprets 2π as the mathematical constant 2π. Secondly, yes, that’s a π! Julia uses allows unicode to be entered (In Juno you enter the latex pi and hit tab), and some of them are set to their appropriate mathematical constants. This is only cosmetic, but in the long-run I can see this as really beneficial for checking the implementation of equations.

However, not everything is rosy in Julia land. For one, many packages, including the FEM package I am working with, are in MATLAB. Luckily, the interfacing via MATLAB.jl tends to work really well. In the first post I showed how to do simple function calls. This did not work for what I needed to do since these function calls don’t know how to pass a function handle. However, digging into MATLAB.jl’s package I found out that I could do the following:

  put_variable(get_default_msession(),:node,node)
  put_variable(get_default_msession(),:elem,elem)
  put_variable(get_default_msession(),:u,u)
  eval_string(get_default_msession(),"sol = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))/(8*pi*pi);")
  eval_string(get_default_msession(),"Du = @(x)([cos(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))./(4*pi) -sin(2*pi.*x(:,1)).*sin(2*pi.*x(:,2))./(4*pi)]);")
 
  eval_string(get_default_msession(),"h1 = getH1error(node,elem,Du,u);");
  eval_string(get_default_msession(),"l2 = getL2error(node,elem,sol,u);");
  h1 = jscalar(get_mvariable(get_default_msession(),:h1));
  l2 = jscalar(get_mvariable(get_default_msession(),:l2));

Here I send the variables node, elem, and u to MATLAB. Then I directly evaluate strings in MATLAB to make function handles. With all of the variables appropriately in MATLAB, I call the function getH1error and save its value (in MATLAB). I then use the get_mvariable to bring the result into MATLAB. That value is a value of MATLAB type, and so I use the MATLAB.jl’s conversion function jscalar to then get the scalar result h1. As you can see, using MATLAB.jl in this fashion is general enough to do any of your linking needs.

For very specialized packages, this is good. For testing the ported code for correctness, this is great. However, I hope not to do this in general. Sadly, every once in awhile I run into a missing function. In this example, I needed accumarray. It seems I am not the only MATLAB exile as once again Julia implementations are readily available. The lead Julia developer Stefan has a general answer:

function accumarray2(subs, val, fun=sum, fillval=0; sz=maximum(subs,1), issparse=false)
   counts = Dict()
   for i = 1:size(subs,1)
        counts[subs[i,:]]=[get(counts,subs[i,:],[]);val[i...]]
   end 
   A = fillval*ones(sz...) 
   for j = keys(counts)
        A[j...] = fun(counts[j])
   end
   issparse ? sparse(A) : A
end
  0.496260 seconds (2.94 M allocations: 123.156 MB, 8.01% gc time)
  0.536521 seconds (2.94 M allocations: 123.156 MB, 8.83% gc time)
  0.527007 seconds (2.94 M allocations: 123.156 MB, 9.41% gc time)
  0.544096 seconds (2.94 M allocations: 123.156 MB, 9.76% gc time)
  0.526110 seconds (2.94 M allocations: 123.156 MB, 12.22% gc time)

whereas Tim answer has less options but achieves better performance:

function accumarray(subs, val, sz=(maximum(subs),)) 
    A = zeros(eltype(val), sz...) 
    for i = 1:length(val) 
        @inbounds A[subs[i]] += val[i] 
    end 
    A 
end

Timings

  0.000355 seconds (10 allocations: 548.813 KB)
  0.000256 seconds (10 allocations: 548.813 KB)
  0.000556 seconds (10 allocations: 548.813 KB)
  0.000529 seconds (10 allocations: 548.813 KB)
  0.000536 seconds (10 allocations: 548.813 KB)
  0.000379 seconds (10 allocations: 548.813 KB)

Why is Julia “Missing” So Many Functions? And How Do We Fix It?

This is not the first MATLAB function I found myself missing. Off the top of my head I know I had to get/make versions of meshgrid, dot(var,dimension) just in the last week. To the dismay of many MATLAB exiles, many of the Julia developers are against “cluttering the base” with these types of functions. While it is easy to implement these routines yourself, many of these routines are simple and repeated by mathematical programmers around the world. By setting a standard name and implementation to the function, it helps code-reusability and interpretability.

However, the developers do make a good point that there is no reason for these functions to be in the Base. Julia’s Base is for functions that are required for general use and should be kept small in order to make it easier for the developers to focus on the core functionality and limit the resources required for a standard Julia install. This will increase the number of places where Julia could be used/adopted, and will help ensure the namespace isn’t too full (i.e. you’re not stepping on too many pre-made functions).

But mathematicians need these functions. That is why I will be starting a Julia Extended Mathematical Package. This package will be to hold the functions that are not essential language functions, but “essential math language” functions like you’d find in the base of MATLAB, R, numpy/scipy, etc., or even just really useful routines for mathematical programming. I plan on cleaning up the MATLAB implementations I have found/made ASAP and putting this package up on github for others to contribute to. My hope is to have a pretty strong package that contains the helper functions you’d expect to have in a mathematical programming language. Then just by typing using ExtendedMath, you will have access to all the special mathematical functionality you’re used to.

Conclusion

As of now I have a working FEM solver in Julia for Poisson’s equation with mixed Neumann and Dirichlet boundary conditions. This code has been tested for convergence and accuracy and is successful. However, this code has some calls to MATLAB. I hope to clean this up and after this portion is autonomous, I will open up the repository. The next stage will be to add more solvers: more equations, adaptive solvers, etc., as I go through the course. As, as mentioned before, I will be refactoring out the standard mathematical routines and putting that to a different library which I hope to get running ASAP for others to start contributing to. Stay tuned!

The post Julia iFEM3: Solving the Poisson Equation appeared first on Stochastic Lifestyle.

Julia iFEM 2: Optimizing Stiffness Matrix Assembly

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem2/

This is the second post looking at building a finite element method solver in Julia. The first post was about mesh generation and language bindings. In this post we are going to focus on performance. We start with the command from the previous post:

node,elem = squaremesh([0 1 0 1],.01)

which generates an array elem where each row holds the reference indices to the 3 points which form a triangle (element). The actual locations of these points are in the array node, and so node(1) gives the points in the (x,y)-plane for the ith point. What the call is saying is that these are generated for the unit square with mesh-size .01, meaning we have 10201 triangles.

The approach to building the stiffness matrix for the Poisson equation is described here. The general idea is that span our vector space by a basis of hat functions phi_{i}, and the so the stiffness matrix is found by the inner product (integral) between these basis functions. This translates to solving for the area of the triangles where two hat functions overlap, which we can do exactly since we chose the basis to be sufficiently nice. However, since the hat functions are zero except in a small range, most of these inner products are zero, meaning the resulting stiffness matrix is sparse. Our goal is to produce this sparse matrix as efficiently as possible. Lets get to it!

Building the Julia Version: Local Stiffness

The first function we need is titled localStiffness, which evaluates the inner product to give the “stiffness for one triangle”. In MATLAB the code is of the form:

function [At,area] = localstiffness(p)
    At = zeros(3,3);
    B = [p(1,:)-p(3,:); p(2,:)-p(3,:)];
    G = [[1,0]',[0,1]',[-1,-1]'];
    area = 0.5*abs(det(B));
    for i = 1:3
        for j = 1:3
            At(i,j) = area*((BG(:,i))'*(BG(:,j)));
        end
    end
end

It takes in a vector p of three points, solves for the area at the reference triangle, and transforms the area appropriately to give the stiffness for the triangle defined by the points of p. I want to make a few points about this code. First of all, it employs a “trick” for solving for the dot product. That is, it uses the transposed vector times another vector. I quote trick because to a mathematician, this is simply the definition of the dot product, and so it only seems natural to use it like this in MATLAB. However, two things to point out. First of all, in Julia, such an operation does not return a scalar, but a one dimensional vector. This in some cases will give unexpected errors, so it should be avoided. Not only that, but it is also inefficient. In both MATLAB and Julia, matrices are stored column-wise, that is they are stored as a array of pointers where the pointers go to an array of columns. Thus to access a row matrix, both MATLAB and Julia would have to access the pointer and then go to the array at which it points (a size 1 array), and take the value there. Notice this is an extra step. Therefore it’s more efficient to keep the vectors column-wise. (In reality, the vectors here are so small that it won’t make a difference, but this justifies that the change we will do for the Julia code in a performance-sense).

Two other small changes were made to this code. For one, Julia throws an error at the definition of G since it cannot read the transpose calls inside of an array declaration. This is easily fixed by simply transposing after the creation instead. Lastly, we change the pre-allocation of At from zeros to Julia’s general constructor. This is slightly more performant since it will allocate the space without doing an initial re-write step, saving the time it would take to loop through and set each value to zero. The result is the following:

function localstiffness(p)
  At = Array{Float64}(3,3);
  B = [p[1,:]-p[3,:]; p[2,:]-p[3,:]];
  G = [1 0 ; 0 1 ; -1 -1]';
  area = 0.5*abs(det(B));
  for i = 1:3
    for j = 1:3
      At[i,j] = area*dot(BG[:,i],BG[:,j]); 
    end
  end
  return(At,area)
end

which is ever so slightly more efficient, but in reality the same and just tweaked for the quirks.

Building the Julia Version: Matrix Assembly

Now we need to loop through each triangle and sum up the inner product between each pair of basis functions over each triangle. The intuitive code is:

function assemblingstandard(node,elem)
  N=size(node,1); NT=size(elem,1);
  A=zeros(N,N);
  for t=1:NT
    At,=localstiffness(node[vec(elem[t,:]),:]);
    for i=1:3
      @simd for j=1:3
        @inbounds A[elem[t,i],elem[t,j]]=A[elem[t,i],elem[t,j]]+At[i,j];
      end
    end
  end
  return(A)
end

I discussed previously the use of macros to speed up code without changing its style. The main problem that had to addressed in porting the code to Julia was that Julia will only take a vector for array referencing when the colon operator is used. Therefore since elem[t,:] returns a 1×3 Matrix of indices for the points associated with triangle t, once again a 1×3 Matrix is not an array in Julia so it throws an error. This is easy to fix by wrapping it in the vec() command, which others have tested to be the fastest method for conversion, and will actually do some fanciness in the background in order to not have to copy the array. This means that the cost of using vec is so small that I will use it liberally as a fix in these cases. Notice that within the loop vec is not required to reference A. This is because the issue only occurs when the colon operator is present.

However, since At is usually zero, we can improve this a lot by instead generating vectors to build a sparse matrix. What we will instead do is save the (i,j) pairs where the value should be stored, and the value, and use the sparse command to reduce. The sparse command will automatically add together the values from repeated (i,j) indices, effectively performing the update we had before. This gives us the code:

function assemblingsparse(node,elem)
  N = size(node,1); NT = size(elem,1);
  i = Array{Int64}(NT,3,3); j = Array{Int64}(NT,3,3); s = zeros(NT,3,3);
  index = 0;
  for t = 1:NT
    At, = localstiffness(node[vec(elem[t,:]),:]);
    for ti = 1:3, tj = 1:3
        i[t,ti,tj] = elem[t,ti];
        j[t,ti,tj] = elem[t,tj];
        s[t,ti,tj] = At[ti,tj];
    end
  end
  return(sparse(vec(i), vec(j), vec(s)));
end

Notice here that vec() needed to be used to build the sparse matrix because the easiest way to index within the loop was to use a 3-dimensional array and then flatten it via vec(). Notice a neat Julia trick where you can define multiple for loops in one line: “for ti = 1:3, tj = 1:3” is two loops and works as you’d expect. With some extra work we can get rid of the loop over the triangles by performing the calculations from localstiffness vector-wise, which gives us the vectorized form:

function assembling(node,elem)
  N = size(node,1); NT = size(elem,1);
  ii = Vector{Int64}(9*NT); jj = Vector{Int64}(9*NT); sA = Vector{Float64}(9*NT);
  ve = Array{Float64}(NT,2,3)
  ve[:,:,3] = node[vec(elem[:,2]),:]-node[vec(elem[:,1]),:];
  ve[:,:,1] = node[vec(elem[:,3]),:]-node[vec(elem[:,2]),:];
  ve[:,:,2] = node[vec(elem[:,1]),:]-node[vec(elem[:,3]),:];
  area = 0.5*abs(-ve[:,1,3].*ve[:,2,2]+ve[:,2,3].*ve[:,1,2]);
  index = 0;
  for i = 1:3, j = 1:3
     @inbounds begin
     ii[index+1:index+NT] = elem[:,i];
     jj[index+1:index+NT] = elem[:,j];
     sA[index+1:index+NT] = sum(ve[:,:,i].* ve[:,:,j],2) ./(4*area); # Replacing dot(ve[:,:,i],ve[:,:,j],2)
     index = index + NT;
     end
   end
  return(sparse(ii,jj,sA));
end

Here I note that Julia’s dot product will not act on matrices, only vectors. In order to do the row-wise dot product like we would do in MATLAB, we can simply use .* and sum the results in each row.

Dispelling a Myth: Vectorized Julia Rocks!

The most common complaint about Julia that people tend to have is that, in many cases, the code which gets the most performance is the de-vectorized code. “But the vectorized code can be so beautiful! Why would I want to change that?”. This myth seems to come from the alpha days or really bad tests, but it doesn’t seem to die. Instead, what I wish to show here is that vectorized code is also faster in Julia. First we run some basic timings:

@time assemblingstandard(node,elem);
@time assemblingsparse(node,elem);
@time assembling(node,elem);
 
0.538606 seconds (2.52 M allocations: 924.214 MB, 12.79% gc time)
0.322044 seconds (2.52 M allocations: 137.714 MB, 21.53% gc time)
0.015182 seconds (775 allocations: 28.650 MB, 16.97% gc time)

These timings pretty stability show this pattern. All of the methods were tried with parallelization and simd options with either no speedup or it being detrimental given the problem size. What this shows is that, out of the intuitive forms for solving these equations, the vectorized form was by far the fastest. The reason is that this is a highly vectorizable problem, whereas I discussed before the limitations that can cause vectorization to lose to devectorization.

But how does this fare against MATLAB? The “same” code was run in MATLAB (this is from iFEM, an optimized library Professor Long Chen) which gives the results:

tic; assemblingstandard(node,elem); toc;
tic; assemblingsparse(node,elem); toc;
tic; assembling(node,elem); toc;
 
Elapsed time is 0.874919 seconds.
Elapsed time is 0.698763 seconds.
Elapsed time is 0.020453 seconds.

To ensure the time difference between the vectorized versions, we had both problems solve it in a loop:

@time for i = 1:1000
  assembling(node,elem);
end
 
 9.312876 seconds (821.25 k allocations: 27.980 GB, 21.32% gc time)

vs MATLAB:

   tic; 
   for i=1:1000 assembling(node,elem); 
   end
   toc;
 
Elapsed time is 19.221982 seconds.

Thus, in line with what this coder found with R, the vectorized code ran more than twice as fast in Julia.

The Julia Optimization Mentality

This was a fun little exercise because I had no idea how it would turn out. Quite frankly, when I started I thought MATLAB would slightly edge out Julia when running such vectorized code. However, Julia continues to impress me. The only major problem that I had this time around was finding out to use the vec() function to deal with “1-dimensional matrices”, but once that was found I was able to get Julia to be faster than MATLAB, even though I know much more about MATLAB and this package itself is quite well optimized.

I think I should end on why I find the Julia philosophy compelling for scientific computing. The idea is not that “you have to devectorize to get the best code”, though there are situations where doing so can dramatically increase your speed. The idea is that you don’t have to contort yourself to vectorization to make everything work. In MATLAB, R, and Python, you have to vectorize in order to make your code to ever finish. That is not the case in Julia. Here, just write the code that seems natural and it will do really well. In this case, vectorized code was natural, and as you could see we got something that was even faster. To do better in MATLAB, we would at this stage have to start writing C/MEX code. In Julia, we could expand out the loops, play with adding SIMD calls, caching, etc. directly in the Julia language.

For scientific computing where we just want code that’s good enough to work, you can see it’s easier to get there with Julia. If you need to optimize it to be part of a library, you can optimize and devectorize it within Julia without having to go to C (many times by just adding macros throughout your code). Will it be as fast as C? No, but many tests show that you’ll at least get within a factor of two so, for almost every case, you might as well just code it in Julia and move on. Each blog post I do I am getting more and more converted!

The post Julia iFEM 2: Optimizing Stiffness Matrix Assembly appeared first on Stochastic Lifestyle.