I'm the primary author of many Julia libraries geared toward "architecture-specific optimizations", including LoopVectorization.jl and Octavian.jl.
For BLAS-like operations, one of the most important optimizations LoopVectorization.jl does is "register tiling".
While CPUs may have a huge number of actual registers (used for register ...
Some things I can think of:
use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY
use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system.
you are solving with ode15s, which is an implicit ...
Abdullah, thank you for the plug for our materials. We have repackaged these as a Massive Open Online Course (MOOC) on edX titled "LAFF-On Programming for High Performance". It is free for auditors. For info, see http://ulaff.net.
The premier open source implementation of the BLAS is the BLAS-like Library Instantiation Software (BLIS). There is ...
The simple example as posted wasn't stiff, but I put it together in Julia anyways for show. I modified it to be a system of 2 PDEs with N=1200 and get in the 10's of ms
using ModelingToolkit, LinearAlgebra, BenchmarkTools, DifferentialEquations
# Setup matrices
N = 1200
mat1 = hcat(zeros(N),Tridiagonal(ones(N-1),-2*ones(N),ones(N-1)),zeros(N))
mat1[1,1] = 1
As mentioned, netlib BLAS is not at all optimized, but it is definetly the "refblas". Using IKML, ACML, OpenBLAS or "your vendor" BLAS, you are (somehow) assured, that the results of the operation of the optimized BLAS is equal to the "refblas" up to a known error. Take into care that: vendors (intel, amd, nvidia, ...) try hard ...