What is fast way do a tensor contraction for two indexes by BLAS(dgemm)?

Question

I want to calculate tensor contraction. X[a,b,e] = Y[a,b,c,d] * Z[c,d,e] Here, I want to contract two indexes c and d together. My idea is:

1. reshape Y=[a,b,c,d] ==> Y[a,b,m] , m=c*d; Z[c,d,e] ==> Z[m,e]
2. X = dgmm(Y,Z)

Reshaping costs some time. Is there any other more efficient way to contract two indexes by fortran?