No example sorry, but at some point I was working on a model with a sparse transition matrix (mostly diagonal) and found good speedups by writing a wrapper around the matrix_exp function to check for zeroes in the relevant off diagonals and compute accordingly. I just ditched it from my own code as part of ‘chuck everything overboard it won’t compile on 32 bit and it’s cran submission time’, but here it is ;)

```
int[] checkoffdiagzero(matrix M){
int z[rows(M)];
for(i in 1:rows(M)){
z[i] = 0;
for(j in 1:cols(M)){
if(i!=j){
if(M[i,j] != 0){
z[i] = 1;
break;
}
}
}
if(z[i]==0){ //check rows
for(j in 1:rows(M)){
if(i!=j){
if(M[j,i] != 0){
z[i] = 1;
break;
}
}
}
}
}
return z;
}
matrix matrix_exp(matrix M){
matrix[rows(M),rows(M)] out;
int z[rows(out)] = checkoffdiagzero(M);
int z1[sum(z)];
int z0[rows(M)-sum(z)];
int cz1 = 1;
int cz0 = 1;
for(i in 1:rows(M)){
if(z[i] == 1){
z1[cz1] = i;
cz1 += 1;
}
if(z[i] == 0){
z0[cz0] = i;
cz0 += 1;
}
}
if(size(z1) > 0) out[z1,z1] = matrix_exp(M[z1,z1]);
out[z0,] = rep_matrix(0,size(z0),rows(M));
out[,z0] = rep_matrix(0,rows(M),size(z0));
for(i in 1:size(z0)) out[z0[i],z0[i]] = exp(M[z0[i],z0[i]]);
return out;
}
```