A 5000 by 5000 matrix contains 25 million elements. That is 200 million bytes of 64 bit doubles, together with 5001 perl array header structures. You need three of those monsters to do an array multiplication. Do you suppose you're hitting swap?

Since most N-element matrix multiplication schemes involve O(N**3) multiplications, you're on the bad side of a scaling problem, as well.

Your strategy of serializing to disk looks pretty good in light of that. You may want to transpose the right hand matrix in order to get the columns as a single array ref.

The best cure would be less brute force. Try to reduce your problem to a block diagonal form if you can. Think about transformations that may simplify things. Would approximate knowledge of the eigenvectors help?

In any case, PDL is superior to pure perl math packages for speed and flexibility.

After Compline,
Zaxo