Currently the implementation of block diagonal function operations is similar to the Matlab version using module of functions. It will be better to reimplement all related functions using the class method. It is also more convenient to be able to take a vector as input in functions such as blk_diag_apply and blk_diag_solve.