Hi @gbarter, thanks for this nice code! By skimming through it I noticed that you have computed a couple of matrix inverses. Generally, this has bad numerical properties, as it's both slow and inaccurate. You should instead solve a linear system.
I suggest to:
At raft_model L399 replace np.matmul(np.linalg.inv(M_tot), C_tot) with np.linalg.solve(M_tot, C_tot).
At raft_model L540 replace np.matmul(np.linalg.inv(Z[:,:,ii]), F_tot[:,ii] ) with np.linalg.solve(Z[:,:,ii], F_tot[:,ii]).
Furthermore, if M_tot and Z are reasonably sparse, you might get an even faster solution using spsolve.
By the way, the routines from scipy.linalg might be faster than the ones from numpy on some architectures.
Hi @gbarter, thanks for this nice code! By skimming through it I noticed that you have computed a couple of matrix inverses. Generally, this has bad numerical properties, as it's both slow and inaccurate. You should instead solve a linear system.
I suggest to:
np.matmul(np.linalg.inv(M_tot), C_tot)
withnp.linalg.solve(M_tot, C_tot)
.np.matmul(np.linalg.inv(Z[:,:,ii]), F_tot[:,ii] )
withnp.linalg.solve(Z[:,:,ii], F_tot[:,ii])
.Furthermore, if
M_tot
andZ
are reasonably sparse, you might get an even faster solution usingspsolve
.By the way, the routines from scipy.linalg might be faster than the ones from numpy on some architectures.